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ABSTRACT 


The aim of the paper is to analyze data obtained during sounding the Gdansk Bay sea bed by means of the 
parametric echo-sounder. The accuracy of the sea bottom structure investigation needs correct configuration 
of research equipment and proper calibration of peripheral devices (GPS, heading sensor, MRU-Z motion 
sensor and navigation instruments which provide necessary data to bathymetrical measurement system, 
enabling its work with full capability. The parametric sub-bottom profiler delivers two types of data which 
include different information: envelope and pure echo signal sampled without processing. The first is used 
to present echograms in real time and the second one is stored during sounding and the obtained data 
can be analyzed by using post-processing software. Differences in the results are presented and discussed. 
And, are described also criteria taken into account during determination of the measurement areas and 
selection of different configurations of the software according to the information to be obtained from 
trials. Complementary information can be achieved from processing the data by means of the procedures 
offered by Matlab software, which allows for a clear 3D presentation of the results. The consideration are 
illustrated by echograms taken in the natural conditions. 


Keywords: sea bottom acoustic; sub-bottom profiler; sea bed structure 


INTRODUCTION 


A bathymetrical measurement system fitted with SES-2000 
sub-bottom profiler was installed on a small survey vessel of 
10.5 m in length. Parametric echo-sounder antenna was installed 
on starboard and EM3002 multibeam transducer - on port, 100 cm 
below water surface on special mounting arms. Additional 
navigation devices were tested in different localizations and 
finally mounted in the places best suited to their functions. 
GPS was installed in the vessel’s centre line, close to deck to 
minimize speed and position errors. MRU-Z motion sensor 
was fitted near the vessel’s gravity centre. Application of the 
SES-2000 Standard - parametric sub-bottom profiler makes 
it possible to observe structure of the upper layers of seabed. 
Parametric sources used in studying the seabed allow - as 
a result of the minimum impact of two acoustic beams radiated 
into the water in the same area - to generate low frequency 
waves (of a few kHz only) of the geometric properties similar to 
the characteristics of the so-called primary waves. This allows 
for observation of stratification of the seabed with relatively 
high accuracy. In addition, the use of the equipment intended 
for the so-called electronic stabilization, in this particular case 
the MRU-Z system, allows to eliminate the sounding beam 


movement associated with the movement of the antenna. After 
mounting the sounding devices, measurements were done to 
define lever arms and compensation parameters for each of 
the bathymetrical units. Signals from GPS and heading and 
motion sensor are distributed through multiplying modules 
to different devices. Trials on Gdansk Bay with calibrated 
instruments on small measurement vessel were made and 
interesting solutions connected with processing and sounding 
were quoted. At the beginning of this article, is presented and 
discussed information obtained during trials on the water 
tank 30 m long, 3 m wide and 1.5 m deep, owned by Gdansk 
University of Technology. 


METHOD OF DATA PROCESSING 


The analysis of data obtained during trials is preceded by 
information on characteristic parameters of antenna (directivity, 
sounding pulses) and methods of data processing prior to 
echogram visualization. Sounding pulse (of 4kHz frequency 
and two-pulse length) measured at 4 m distance from transducer 
(the near measurement point II marked in Fig. 2), is plotted in 
Fig. 1. The calculated depth resolution of such configuration 
amounts to 37.5 cm. 
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Fig. 1. Sounding pulse measured at 4 m distance from transducer 


Measurement Point 11 Measurement Point I — 


Fig. 2. Schematic view of the water tank 


The data visualization method are presented on the example 
of sounding the concrete wall. Plots shown in the figures below 
were recorded with the use of Tektronix oscilloscope connected 
to analogue output of the echo-sounder, which provides a filtered 
and amplified signal just before sending it to analogue-digital 
converter (such signal is shown in Fig. 3). Fig. 4 illustrates 
the envelope calculated from pure echo signal. 


Amplitude [mV] 


Time [ms] 


Fig. 3. Echo of sounding pulse 


The measured signal length is 1.5 m, but not all its fragments 
should take part in visualization, that is shown in Fig. 5 
where the gray flat area is the threshold, i.e. LF Min Level 
for small echoes and noise. Changing this value can strongly 
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Amplitude |A/Amax| 


Time [ms] 


Fig. 4. Envelope of received pulse 


influence visualization model and help examination of received 
information. The typical echogram which describes result of 
the sounding with all cut-off filters set to minimum values, is 
shown in Fig. 6. 


Amplitude |A/Amax| 
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Fig. 5. 3D visualization performed by using the Matlab software 
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Fig. 6. 2D echogram received from the Matlab software 


The data received from the stationary measurements 
in 30 m water tank, were used during analyzing the results 
from the trials at Gdansk Bay. It was important to set depth 
resolution for different sounding pulses depending on their 
frequency (4, 6, 8, 10, 12 - 15 kHz) and number of pulses 
which represent number of cycles. Increasing number of pulses 
increases energy of transmitted signal (better signal to noise 
ratio) but decreases the depth resolution: 


Depth_resolution = t*c/2 


Where: 
t = 
1G => 


time, 
sound velocity in water 


Fig. 7. Measurement layers with SES data file possible 
to be processed after recording 
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Fig. 8. Echo and envelope of received signal 


The result obtained from sounding two artificial sediment 
layers is shown in Fig. 8. The two artificial sediment layers 
were arranged in the configuration shown in Fig. 7, where their 
mutual distance was 35 cm. Looking at the result, one can notice 
a characteristic point in which a jump in the envelope (shown 
on echogram as free space) is formed during calculating the 
envelope. The free space length depends on threshold level and 
type of visualization method, but the length is correlated with 
free space between frames. 


DATA RECEIVED FROM SOUNDING 


Comparison of pure signal 
and lower resolution envelope 


Data from the sounding at Gdansk Bay where penetration 
depth of sediments reached 40 meters (the starting point is 
set at 50 m below sea surface) are presented in Fig. 9. 4 kHz 
frequency and 2 pulses length signal is used. The two following 
charts describe comparison of the data: 

- SES data — contain 480 samples used to draw echograms 
in real time. 

- RAW data — contain 5120 samples which are stored in 
parallel. 


-B= 


o ea T 


-a--G 


Fig. 9. Comparison of *.SES and *.RAW data 
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The interesting comparison of the data was expanded in 
the upper part of Fig. 9. Difference in their quality can be 
observed, which clearly shows that the slow sampling brings 
errors connected with resolution, causes losing information 
about e.g. rapid drops in envelope, that can be very important 
for searching the objects buried not deeply. The envelope 
calculated from the recorded pure echo signal is shown in the 
lower part of Fig. 9 — its expanded part shows all details without 
losing information. After using proper post-processing tools 
and *.raw data files small changes in signals can be observed 
during analyzing full echo runs. 


Filtering pure signal and data presentation 
on echogram 


The below described results are obtained from Matlab 
software where appropriate algorithms for data analyzing are 
implemented. Additionally, interesting results of sounding were 
selected to show rapid change in sediment properties. The first 
370 m recorded on echogram show the sea bottom located in 
depth of 54 m and sediment layers up to 85 m. After sounding 
380 m distance we can observe fast change of penetration 
depth which decreases to 2.5 m. Looking forward the nearest 
1000 m of the visualized echogram one can observe a stronger 
echo reflected from the bottom and no information from the 
sub-bottom zone. 

3D visualization without filtering 


— 


Recon ATEN 


Distance [m] 


Fig. 10. Raw data displayed in the form of echogram - pure signal 


70 
Depth [m] 


Having known the sounding pulse and sampling frequency 
it was possible to filter signals by using the band-pass filter. 
The result of the filtering is shown in Fig. 11, where is possible 
to observe the clear echogram with more differentiated echo 
levels in the three marked zones. 


3D visualization with filtering 


Distance [m] 
Fig. 11. Raw data displayed in the form of echogram - signal passed 
through band-pass filter 


The time series of the reflected pulses received after 
transmitting sounding pulse were plotted on three charts (Fig. 
12). The first one visualizes results from the sounding zone I. 
The analyzing of the signal did not give clear classification 
of echo reflected from layers - signals are continuous, but 
after forming the envelope characteristic hooks (which arise 
during impact of two overlapping reflected sounding pulses) 
allow to plot echograms with visible boundary between layers. 
The signal from the middle part of bottom sounding (zone IT) 
exhibits echo from the first layer located on a higher level than 
in the zone I, that can be interpreted for instance as another type 
of bottom material imposed on a muddy bottom. 

Looking at the third time series one can observe that in this 
case the first and second echo show much higher level than in 
the two first zones. To get proper interpretation of sounding 
results (especially in the third case - possible simulation of such 
layer arrangement) a correlation with the data received from 
trials on the 30 m water tank is necessary where differential 
sounding pulses were a.o. measured and the sounding of two 
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Fig. 12. Time series of raw signals 
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layers side by side with different types of bottom sediments 
placed in special frames, was performed. 

4 kHz sounding pulse is shown in Fig. 13. The sounding 
pulse is nearly the same as the reflected signal from zone HI, 
that can be checked during analyzing Fig. 12. The first echo 
is a little fuzzy, that can be interpreted as reflection from 
a muddy/clay material but the second one corresponds directly 
to the transmitted signal. 

Sounding pulse 4kHz 1 pulse 


trials in the water tank (lower figure) 


To compare echo power level, respective soundings from 
each zone, after filtration, were averaged and displayed in 
Fig. 14 in the form of artificial echogram. 

3D visualization arfter processing - with filter 


70 
Depth [m] 
80 


500 
Distance [m] 


90 9 
Fig. 14. Comparison of echo strength 


The echogram confirms that the deepest bottom penetration 
has lower echo level along whole path. Strength of the signal 
reflected from the area near zone I is greater but not as great as 
the signal coming from zone III. During sea trials several points 
rapid change of sediment penetration depth were localized, that 
gives opportunity to mark the points on map and apportionment 
border on bay (red colour line) whose shape is regular and can 
confirm real natural properties of the sediments in this area. 

The deepest penetration zone is marked yellow and the red 
area corresponds to place where investigation of sediments 
is reduced to that of the high power echo reflected from the 
first layer (the sediment samples, taken not from the whole 
coastal zone, were analyzed during water tank laboratory 
measurements, where they were placed in the frames to be 
sounded with low frequency pulses emitted by the same unit 
as during sea bottom trials). 


df | 


Ri 


— N å 
Area with hard bottom - 
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Fig. 15. Red line - border of penetration 


After expanding an interesting fragment of the echogram 
(Fig. 16) it is possible to calculate a distance from which the 
penetration depth decreases from 20 to 5 m. The distance 
is near 13 m which corresponds to a fast change in bottom 
properties. However such geological case can also appear 
under natural conditions and this change does not result from 
one sounding ping to another so it is not due to a change in 
sounding parameters, like e.g. frequency. 


Fig. 16. Expansion of rapid change of penetration depth 


The sounding results which show rapid change of 
penetration depth are correlated with results obtained by 
another researchers and published on relevant websites. The 
results obtained by these authors are equivalent to echograms 
produced by different echo-sounders, like e.g. the TOPAS 40 
- results titled “Results from shallow gas area outside Namibia”, 
published on Kongsberg website. 
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Fig. 17. Kongsberg’s TOPAS 40 echogram 


It is necessary to point out that during measurements 
a compensation of roll (beam stabilization) and heave motions 
was realized with the use of the MRU-Z device. But it is 
important to remember of influence of not compensated pitch 
motion, which depends on sea state and can be considerable. It 
is easy to confirm the influence during analyzing the echogram 
shown in Fig. 18 where bottom sediments are the same but pitch 
motion makes sounding pulses hitting bottom at different angles, 
that changes echo amplitude. These results were taken during 
high surge motion in order to get an example of such data. The 
typical sounding operation was realized during calm weather. 


Fig. 18. Example of different levels of echo reflected from the same sediment 


METHODS OF POST-PROCESSING 
VISUALIZATION 


During searching for difference between sediment properties 
or checking a mud or silt layer covering sandy bottom it is 
important to set, in a proper way, threshold for minimum values 
to be visualized on echogram, in order to make substantial 
differences visible. Additionally, the algorithm which defines 
distribution of colours should be set to logarithmic or linear 
mode. By using the Matlab as programming software to 
visualize data, two examples taken during sounding two areas 
close to each other, are shown below. 

Each echogram contains two parts — on the left side are 
visualized sounding results where penetration depth was up 
to 3 m, and on the right side - echoes reflected from the layers 
located 15 m under the bottom. The colour scale was set to 
linear mode. In Fig. 19 the threshold was set to minimum to 
make all details visible, but usually it is difficult (especially 
in logarithmic colour scale mode) to distinguish a mud layer 
placed on hard bottom. The threshold level was set to the value 
which cuts off envelope picks resulting from low power echoes, 
shown in Fig. 20. It can be observed that the three thick layers 
which were first the bottom layer, now are covered with 3 m 
thick sediment. The same comparison can be seen in 3D view 
of envelopes, as shown in Fig. 21. The envelope used in the 
comparison is averaged from sounding records, that makes its 
peaks smooth. Additionally, when data include much noise it 
is possible to use a stacking algorithm (adding consecutive 
echo signals and receiving new, processed data) and smoothing 
function which decreases influence of noise on the results. 
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Depth 50 m 


Depth 50 m 


Fig. 20. Echogram with high threshold value 


0m depth 


5 


Fig. 21. 3D presentation of data by using Matlab software 


Operator should select proper configuration of the software 
intended for co-operating with parametric echo-sounder, 
depending on the type of research to be performed, namely: 
a) examination of sediments layers, 

b) searching for a buried object. 


Having known that the envelope contains constant number 
of samples, not depending on range, when searching for a small 
buried object located at constant depth on flat bottom, the start 
sampling range can be set 1 +2 m above the bottom and 5 m in 
its value, i.e. the smallest one. It will give the highest resolution 
which provides the sampling step of 0.01 m, that makes it possible 
to observe small changes in reflected sounding pulse. When the 
main aim is to investigate the sediment layers up to 30 m below 
the bottom surface, operator should set up the start sampling 
range to be 4+ 5 m above the bottom surface and the sounding 
range should be set to 35 + 40 m. The resulting resolution is not 
very high but satisfactory in the case in question. 

If something interesting is noticed in the seabed structure, 
for example at 10 m depth below the bottom surface, the 
sounding should be repeated with the sampling window set 
up in the way described in the searching for buried object but 
with the start sampling range value adjusted to actual situation 
(it could be set below bottom surface). 


The sounding results obtained from the software intended 
for imaging data in real time are presented in Fig. 22 and 23, and 
the results post-processed with the use of the software dedicated 
to this type of data — in Fig. 24 and 25. In the considered case 
the bottom is located 46 + 48 m below sea surface and the 
penetration range reaches 27 m in the maximum point. 


Fig. 22. Envelope from 480 samples - 6 kHz 2p 
(from SES data - differential view) 


=e 


Fig. 23. Envelope from 480 samples - 6 kHz 2p 
(from SES data - full envelope view) 
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Fig. 24. Envelope from 5120 samples — 6 kHz 2p 
(RAW data - envelope with 15 colours) 
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Fig. 25. Envelope from 5120 samples — 6 kHz 2p 
(RAW data - envelope with 30 gray-scale colours) 


The software dedicated to the sub-bottom profiler offers 
various configurations. Most important is to set up frequencies 
of differential sounding pulses, which can be chosen in the 
range of 4 kHz +15 kHz. The primary frequencies are set up 
near 100 kHz, depending on a differential sounding pulse. 
During work with geological aspects where it is important 
to investigate structure of sediments, the lowest frequency is 
required, that can be checked on the below presented echograms 
(where high differences were used for visualization of the same 
sounded bottom). When a new research area is examined it 
is useful to apply multi-frequency option. It allows to sound 
with the use of three different sounding pulses one by one. An 
example of such research is illustrated in Fig. 12, 13 and 14. 
The results are elaborated in the logarithmic scale of colours 
and the TVG software set to 0.6 dB/m. 


Fig. 26. 5 kHz frequency of echosounder souding pulse 


By using data processing algorithms and Matlab software 
to compare and analyze data which contain information about 
buried objects, is possible. There are many possibilities to 
illustrate processed data in 2D, 3D, horizontal and vertical view. 
Having geographical position linked with sounding data it is 
possible to visualize different models of searched objects. Below 
an example of visualization and measurement of three buried 
targets is shown. Additionally, application of scale of colours, 
where user can define step, threshold, maximum and minimum 
values, provides a wide range of visualization capabilities. 
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Fig. 28. 12 kHz Frequency of echosounder souding pulse 
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Fig. 29. Scale of colours of echosounder souding pulse 
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Fig. 30. Set of default colours used in the Matlab 


SUMMARY 


Interpolating colorspace: | RGB x] 


Immediate apply 


During the post-processing of data different algorithms 
should be applied depending on a type of research (examination 
of sediment layers or searching for buried objects). The data 
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Differential frequency echogram 


Depth on 35,5 m 


Fig. 31. Three buried objects - imaged by using the Matlab 


analyzing should additionally improve operator’s knowledge 
useful for new trials, as to the describing of proper configuration 
of a software intended for cooperating with parametric echo- 
sounder, frequency and length (number of cycles) of pulses to 
be set depending on currently received results. Methods of data 
processing in real time and raw data post-processing provide 
many possibilities. It is necessary to have at one’s disposal an 
algorithm which describes a proper way of effective work with 
the data, even during trials. 
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ABSTRACT 


Annular cavitator with water injection is one of the key parts of the long-range supercavitating vehicle 
powered by water ramjet. In this paper, hydrodynamic properties of annular cavitator are studied 
numerically. The standard k ~ € turbulence model is coupled with the Reynolds Averaged Navier-Stokes 
(RANS) equations to model the natural supercavitation process. The multiphase flow is considered as 
a mixture of varying density and modeled by the mass exchange equations. To fully understand this process, 
numerical simulations were performed for different annular cavitators. Computational Fluid Dynamics 
(CFD) results, including the pressure distribution and forces acting on the cavitator surface, mass flow 
and pressure loss of water injection, various supercavity sizes, were obtained and analyzed. The pressure 
distribution on the cavitator surface was significantly changed which resulted in 4 ~ 6% increase of the 
total drag of the vehicle. The results show that the mass flow and velocity of the injection water is mainly 
dependent on the tube size, while the total pressure loss of the water injection is mostly related to the outlet 
pressure. Supercavity generated by annular cavitator is smaller than that of the discal one. Based on the 
correlation analysis of the supercavity size and other factors, it could be concluded that the contraction of 
the cavity size is mainly caused by the diffluent mass flow of the water injection. 


Key words: supercavitation; computation fluid dynamics (CFD); annular cavitator; water injection 


INTRODUCTION 


Cavitation, which happens when a vehicle travels fast 
enough underwater, brings extra noise, surface erosion, and 
other problems. To avoid cavitation, traditional underwater 
vehicles could hardly travel faster than 70 kn. Supercavitation, 
in contrast, which could reduce the drag by more than 90%, 
has been proved to be a revolutionary way to achieve ultrahigh 
speed underwater [1]. Due to the successful applications of the 
Russian torpedo ‘Shkval’ and the US supercavitating projectile 
‘RAMICS’, supercavitating vehicles have received increasing 
research interests in recent years [2-4]. 

With stunt shape and sharp edge, cavitator generates and 
sustains a supercavity to envelope the whole vehicle during high 
speed motion. Both theoretical and numerical studies have been 
carried out to enhance the understanding of various properties 
of cavitators and the supercavitating flow. 

By applying the theorem of the independence of cavity 
section expansion proposed by Logvinovich [5, 6], pioneering 
studies have been conducted by Russian and Ukrainian 
researchers such that they have been widely used for the 
predicting of the cavity shape [7, 8]. L.Sedov studied the jet 
cavitator which ejects water stream into the flow. In this case 
the body will suffer only half drag when the cavitation number 
is near zero [9]. Using computationally supercavitating flow, 


Owis et al considered the compressibility of supercavitating 
flow and improved a numerical method for both single and 
multi-phase flows [10-12]. Drag force for NACA66 foil and 
a flat plate was investigated by Seif and the CFD model was 
successfully tested [13]. Carried out by Kuklinski, experimental 
studies of hydrodynamic properties of different cavitators 
provided test data for the dynamic modeling of cavitators [14]. 
Studied by Ahn, Shafaghat, and Lin, optimal designs of conical 
or spherical cavitators proved that a proper shape could improve 
performance of the cavitator [15-17]. To achieve long-range 
supercavitating motion, the water ramjet is the best propulsion 
device for supercavitating vehicles. In this case, water injection 
from the cavitator is required, which makes conducting research 
on annular cavitators necessary. 

In this paper hydrodynamic properties of the annular 
cavitator are studied in detail. Based on the Reynolds Averaged 
Navier-Stokes (RANS) equations which is coupled with the 
standard k ~ € turbulence model, the supercavitating flow 
is considered as a varying-density single-phase flow. The 
water injection processes with assumed different tube size 
and outlet pressures are studied. Based on the simulation of 
different water injection conditions, total drag of the annular 
cavitators are obtained. The mass flow and the pressure loss 
of the water injection are studied. Related factors which affect 
the supercavity size are investigated by using the correlation 
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analysis. Experimental results shown in this paper can be 
utilized for the designing of supercavitating vehicles propelled 
by water ramjet. 


MODELING AND COMPUTATIONAL 
APPROACHES 


Governing equations 


Based on three dimensional Reynolds Averaged Navier- 
Stokes equations, the supercavitating flow consisted of water 
and vapour is considered as a single phase of the same pressure 
and velocity field and is modelled by the mass exchange 
equations. The steady supercavitation flow is considered 
isothermal, therefore the energy equation is not considered. 

The continuity and momentum equations of the mixture 
phase are given below: 


O(Pm) 
+V 0 (1) 
at 2(PpW) = 
nM) Vg(p,uxu') = 
Ot (2) 


=—Vp + Vg[u,, (Vu + Vu')]+ Png 


The continuity equation of the vapour phase is shown as 
follows: 


OCP a) 
ôt 


where u is the velocity vector of the mixture, g is the 
gravity vector, p, a and u represent density, volume fraction, 
and viscosity, respectively, subscript 1, v and m represent liquid 
phase, vapour phase, and mixture phase separately, m*, m are 
the source terms caused by vaporization and condensation. 
The mixture property, @,,, can be obtained by: 


+ Vg(p, a,u) = m- m (3) 


m? 


Om = ,0, T ,a,, (4) 
where stands for density, viscosity, and so on. 
At last, volume fraction-conservation equation is as 
follows: 


yta, =l (5) 
Natural cavitation model 
The Rayleigh-Plesset equation which provides the basis for 


the rate equation controlling vaporization and condensation, is 
given as follows: 


2 
— 2 
A gee ee (6) 
dt” 2`dt Pi P,Rg 
where 
Rẹ — represents the gas bubble radius, 
p, — saturated vapour pressure, 
p — the pressure in the liquid surrounding the bubble, 
pı  — the liquid density, 
o  — the coefficient of surface tension between the liquid 
and vapour. 


By disregarding the second-order terms and surface tension, 
the rate of vaporization and condensation are shown below: 


m*=C, 30 nui P y 1Py 2 Py — PvP P< p (7) 
Rg 3 PI 
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m=C. 3a, P, 2p p>p, (8) 
R \3 p 
where: 


Onc — the volume fraction of the nucleation sites, 

C., C.— empirical factors which may be different for the rate 
of vaporization and condensation; usually: C, = 50, 
C. = 0.01. 


Turbulence model 


The standard k ~ e turbulence model is used in this study. 
The turbulence kinetic energy, k, and its dissipation rate, ¢,are 
obtained from the following transport equations: 


2 (Pk) + Ve(p,vk) =Ve ure \Vk]+G,-pré (9) 
ô 
—(p,,k) +V = 
A (Pk) + Ve(Ppve) 


E 
=Vg[(p + ZO Ve] +Ë (ChG — Ce Pme) 


where: 

u, = pC,k?/e — the viscosity of turbulence, 

G,=2u,Vv — the kinetic energy of turbulence, 

Op O; — the Prandtl number of k and e, respectively. 
Again, 6,, = 1.44, 6), = 1.92, C, = 0.09, 
6, = 1.0, 6, = 1.3. 


Settings for simulation 


The vehicle model is illustrated in Fig. 1. The length of the 
vehicle is 2 m, while the diameter of the body is 0.2 m. The 
set diameter of the cavitator is 10 cm. To avoid the unstable 
region of the injection hole on the cavitator, the length of the 
tube is set to 1 m. The bottom of the tube is assumed to be the 
outlet surface. The incoming flow velocity is 100 m/s and the 
environmental pressure is 0.2 MPa. 

The generic CFX code was used to investigate the 
liquid flow around the supercavity. Numerical calculations 
were performed over a 180° radial sector of the field with 
symmetrical boundary conditions. The computational mesh 
corresponding to the schematic geometry is shown in Fig. 2. 
In summary, over 500,000 hexahedral elements were formed 
and non-uniformly distributed within the entire computational 
domain. The transport equations were discretized by the finite 
volume approach. The convection terms were approximated 
by a high-order resolution scheme while the diffusion terms 
- by the second-order central difference scheme. Convergence 
was achieved within 2500 iterations when the RMS (root mean 
square) residual dropped below 10-5. 
Annular Water injection 
cavitator tube 


(10) 


Vehicle body 


Fig. 1. Schematic geometry of vehicle model 


Fig. 2. Computational mesh of the annular cavitator 


To consider the effect of the outflow pressure and 
tube diameter, different simulation cases were performed 
accordingly (Tab. 1). The same condition is assumed for case 
3 and 7. The disk cavitator with the same diameter is selected 
to be the case 0 for comparison purpose. 


Tab. 1. Values of tube diameter and outlet pressure for different CFD cases 


RESULTS AND DISCUSSION 
Pressure distributions and drag forces 
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Fig. 3. Pressure distributions on the cavitator surface for case 1 through 4 
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Fig. 4. Pressure distributions on the cavitator surface for case 5 through 8 


Based on the assumption that the static pressure of the 
bottom of the tube to be 4 MPa, the pressure distributions on 
annular cavitators (case 1through 4), are presented in Fig. 3. In 
the radial direction, the pressure increases on the tube section 
and then decreases on the cavitator. As the diameter of the 
injection tube increases, the stagnation ring of the cavitator 
surface moves outward along the radial direction. The pressure 
in the centre is larger than the outlet pressure, which means that 
certain pressure loss happened during the water injection. As 
shown in Fig. 4, although the outlet pressures are different, the 
position of the maximum pressure is almost the same. It means 


that position of stagnation ring is hardly related to the outlet 
pressure, which only depends on the size of the injection tube. 
Furthermore, the pressure distribution on the ring part of the 
cavitator surface is independent on the pressure distribution 
in the tube section. 


Tab. 2. Forces acting on region of cavitator with different tube diameter 


Ema IN] Increment [%] 


i=) 


30608 
31825.4 
32128 
32350 
32366 
32228 
32316 
32350 
32336 


1 
2 
3 
4 
5 
6 
7 
8 


By integrating the pressure in the radial direction the forces 
acting on the cavitator surface are obtained (Tab. 2). As the 
tube diameter increases (case 1 through 4), the cavitator force 
decreases significantly due to contraction of the cavitator surface. 
Meanwhile, the forces acting on the tube section increase faster. 
As a result, the total forces acting on the vehicle increase by 
about 4 ~ 6%. Furthermore, the larger the tube diameter the 
greater increase of the total force. Total forces in the case 5 
through 8 are still larger than the ones in the case 0, but they 
change slightly as the outlet pressure of the tube varies. Such 
variety means that the contribution of pressure in the central 
tube section results in a relatively small contribution to the total 
force. Outlet pressure of the tube shows a lower significance to 
the total force acting onto the supercavitating body. 


Mass flow and pressure loss 


The water injection mass flows and velocities versus tube 
diameters are illustrated in Fig. 5. The results show that the 
mass flow is almost proportional to the tube section area and 
that the velocity of the injecting flow is less sensitive to the tube 
size. As shown in Fig. 6, both the mass flow and water injection 
velocity decrease as the outlet pressure of the tube increases. 
Compared with the results in Fig. 5, the mass flow is rather more 
dependent on the tube size than on the outlet pressure. 

Due to the sudden contraction of the tube section and the 
friction on the tube wall, pressure loss is unavoidable during 
water injection. Since the water is incompressible and the tube 
diameter is constant along the pipe, the total pressure loss is 
equal to the static pressure decrease of the water. The pipe 
static pressure distribution is illustrated in Fig. 7. In all the 
cases significant pressure drop and turbulence can be observed 
near the inlet to the pipe. It means that the pressure distribution 
near the inlet is highly unstable. In case 1, for the smallest 
tube diameter, pressure loss along the pipe is the fastest in 
comparison with the ones for other tube diameters, while case 4 
shows the opposite result. Total pressure loss is larger when the 
injection tube diameter is relatively smaller. In case 1, the total 
loss is about 1 MPa. Regarding the cases 5 through 8, it can be 
concluded that the pressure loss is proportional to the pressure 
difference between the total pressure of flow field and the outlet 
pressure. The smaller outlet pressure the larger pressure loss. In 
case 5, the total pressure loss achieves 1.38 MPa while in case 
8 the loss amounts only to 0.34 MPa. In summary, by applying 
larger tube size and higher outlet pressure, the total pressure 
loss could be reduced dramatically. 
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Fig. 7. Pressure loss of the water injection versus tube length 


Dimension of supercavities in different cases 


According to the CFD algorithm the multi-phase 
supercavitating flow is treated as a mixture of varying density. 
Thus the supercavity boundary is considered as the position 
where the density of flow is p,,=0.5. The supercavities generated 
by different cavitators are presented in Fig. 8. Clear reentrant 
jet flow can be observed in the tail of supercavities. Despite the 
total force increase on cavitator surface, supercavities generated 
by the annular cavitator are smaller than those generated by the 
disc in case 0. We can conclude that the cavity size is dependent 
not only on the drag force. When passing from case | to case 4, 
as the tube diameter increases, the cavity size becomes smaller. 
However, when passing from case 5 to case 8, the cavity size 
grows due to the outlet pressure increase. 
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Fig. 8. Supercavities generated by different cavitators 
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Fig. 9. Coefficients of correlation between cavity size and other possible 
factors 


To underlie the reason of cavity contraction, a correlation 
analysis was performed. The correlation coefficients between 
cavity size and other possible factors are presented in Fig. 9. 
The result shows that the cavity size is positively correlated 
with the cavitator force and mass flow through the tube. The 
force acting on tube section and the tube size is correlated 
negatively. The correlation coefficient between the cavity size 
and water injection mass flow is equal to 0.999, which indicates 
that the supercavity contraction is caused by the water injection 
diffluence in the flow field. 


CONCLUSIONS 


In this paper a numerical research on the natural 
supercavitation problem is performed with the use of a high- 
speed annular cavitator. To model the supercavitating flow, 
the liquid and vapour phase of the flow field is considered as 


a mixture of varying density. The RANS equations coupled 
with the standard k ~ £ turbulence model are applied to CFD 
simulation. Properties of the flow field are obtained for given 
different injection tube size and outlet pressure conditions. 

Pressure distributions, forces acting on the cavitator surface 
and injection water flow properties are analyzed. The obtained 
results show that the stagnation region on annular cavitator 
moves outward along the radial direction, which results in about 
4 ~ 6% increase of the total drag. The tube diameter affects the 
drag increase more than the pipe outlet pressure. The injection 
water mass flow mainly depends on the tube diameter. Both 
the mass flow and injection velocity increase with the tube 
diameter and decrease with the outlet pressure. The design of 
the water injection system will affect not only the amount of 
injected water but also the total pressure loss. The pressure loss 
is inversely proportional to the tube diameter. 

The annular cavitators form smaller supercavities than the 
discal one of the same diameter. According to the correlation 
analysis it is caused by diffluence of the water injection. 

To better understand the performance of a water- ramjet- 
propelled supercavitating vehicle, it is essential to take it 
into account the drag increase and the cavity size contraction 
caused by the annular cavitator. In such case, motion dynamics 
of a supercavitating vehicle may need some modification. 
Furthermore, in terms of mass flow and pressure loss, the 
obtained results would be very helpful in optimal designing 
the water injection system. 
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ABSTRACT 


A hybrid craft is defined as a combination of different lifting surfaces. The lift force can be of aerodynamic 
or hydrodynamic character. This research deals with the conditions of using planing mono-hulls equipped 
with a couple of wings, called planing-wing hybrid craft. A parametric study is conducted by using Vortex 
Lattice Method to determine feasibility of using an aerodynamic wing mounted to planing craft. It has 
been determined that the wing area should be more than twice as much as the wetted surface of the craft, 
to get an effective wing for planing craft. Furthermore, the wing should be so designed as to obtain its 
high aspect ratio. Maximum lift force contribution resulting from different wing configurations considered 
in this study is less than 35 percent of the total hydrodynamic lift force at an arbitrary condition. On this 
basis it may be concluded that the wing may be effectively used only in special conditions, hence this is 
not recommended for general application. 


Key words: hydrodynamics; ground effect; wing; planing; aerodynamics; VLM 


INTRODUCTION 


People are always looking for increasing productivity and 
improving performance of transportation. Serious limitations on 
increasing speed and gradual development of existing marine 
crafts led to a new generation of craft, known as hybrid crafts. 
The objective of this paper is to examine performance of one 
kind of such crafts, namely the planing-wing hybrid craft 
in which hydrodynamic and aerodynamic lift forces mainly 
supports the craft’s weight at service speed. 

Due to higher density and viscosity of water than air, large 
speed ofa body in water causes its frictional and wave-making 
resistance greater than in case of moving in air. One solution to 
reduce this type of resistance is reducing the wetted surface of 
craft. Although lifting surfaces produce a higher lift force and 
reduce craft’s draught but they may cause also some detrimental 
effects such as cavitation of hydrofoils, degradation of craft’s 
performance in waves, reduction of its stability etc. The 
restrictions may put obstacles on increasing the craft speed. 

WIG crafts (Wing- In- Ground- effect crafts) are the crafts 
which make use of aerodynamic forces at service speed. 
Undoubtedly, the use of aerodynamic forces can be considered 
as a worthy contribution to reaching higher speeds in marine 
transportation. The reason is that a significant lift force can 
act onto acraft moving near water surface, by which WIG 
crafts are able to achieve very high speeds. However, the WIG 
crafts are facing some problems such as structural and dynamic 
complexities, construction cost and safety [1] as well as low 
efficiency of air propulsion system. Another disadvantage of 
WIG craft is inability to perform even very simple manoeuvres 
in planing mode. 
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Some types of high speed crafts, so called the hybrid crafts, 
may use hydrodynamic lift of planing hull, hydrodynamic lift 
of hydrofoils and aerodynamic lift of aerofoils. It means that 
in the hybrid craft the capability of various lifting surfaces are 
used to gain more efficiency and performance quality. 

This study is intended for examining usage of aerodynamic 
forces including ground effect and hydrodynamic forces 
simultaneously to increase the efficiency and performance 
quality. A high lift-to-drag ratio of the aerodynamic foils 
can make a fundamental change in power consumption of 
propulsion system and allow reaching higher speeds. 

The main question is what should be the wing dimensions 
so that a significant part of the craft weight to be supported by 
air lift? Is there any advantage from using air lift instead of 
hydrodynamic lift? 

In this context a literature review shows that there are some 
products which use a combination of different kinds of lift 
force, and that none of them have been theoretically discussed. 
In some cases only patent registration was done. In 1976, 
Shipps [2] introduced a type of hybrid craft, known as "tunnel 
hull", used as a race boat. This design in which two planing 
hulls were applied at the end of a ram wing and considered as 
end plates, could promptly demonstrate its better performance 
than planing mono-hulls; and thus a major step was made in 
designing the race crafts. One of the benefits of this design 
was the use of aerodynamic lift force for bearing craft weight. 
In 1978, Ward [3] presented some results of testing the vessel 
KUDU II (a modified form of the vessel KUDU I designed 
in 1976), which could sail at speed up to 78 knot. In 1996, 
Privalov and Kirillovikh [4] designed a hybrid craft called 
TAP (Transport Amphibious Platform) which consisted of 


two catamaran- like bodies, an aircraft fuselage, aerodynamic 
wing and tail. Its motion in close proximity to water surface or 
ground resulted in forming an air cushion between the hulls. 
Among other merits it allows to develop a very high speed (250 
km/h). In 1997, an ekranocat vessel was presented by Doctor 
[5] during WIG's97 Conference, in which an aerodynamic 
superstructure was arranged on catamaran body, resulting in 
drag reduction by 50 percent, according to the calculations 
and measurements. 

Furthermore, several patents for the use of aerodynamic 
elements in high speed crafts have been registered in the US [6], 
and in all of them utilization of ground effect has been indicated 
as the reason for using the wings. In 2005, a research program 
was started at Cranfield University, UK, whose goal was the 
use of aerodynamic forces for high speed craft [7]. 

Balancing the different forces which act on hybrid crafts 
has been a target of some research projects in recent years. 
Collu eż al. [8] and [9] has considered longitudinal stability of 
hybrid crafts in both static and dynamic states. 

The motivation for this study is based on the attractive 
parameters of a boat called Seaphantom, claimed by its 
manufacturer [10]. The Seaphantom is an innovative hybrid 
craft whose original idea formed in 2000 was registered as 
a US patent [6]. In 2004 was built its prototype on which a set 
of tests and trials was performed. The second unit was built 
in 2007 and tested in January 2008. The craft’s designer and 
manufacturer tried to use different lifting elements to achieve 
the maximum possible speed by minimum power. According 
to the manufacturer’s claim, the 300-hp Seaphantom is able 
to move with the speed of 70 knot which is considered an 
excellent feature of the boat. Its wings have a low breadth 
compared to hull breadth, which is other key point to be 
mentioned about the craft. Tab. 1 shows the Seaphantom 
specification. 


Tab. 1. Seaphantom specification [10] 


Quantity 


Length overall 9.27 
Breadth overall 3.818 
Speed(max) 70 
224 


Power 


Engine type Mercury race 


Outboard 


Engine configuration 


Engine speed 9000 


The purpose of this study is to examine hydrodynamic 
advantages of a hybrid craft which uses both the hydrodynamic 
lift force of a planing body and the aerodynamic lift force 
resulting from air wings which are placed close to the water 
level to benefit from aerodynamic lift due to ground effect. 

To achieve the above mentioned objective, a model 
planing-wing hybrid craft on which different air wings are 
assembled, is selected. A schematic image of the hybrid craft 
is shown in Fig. 1. The lift and drag forces produced by main 
body, wings and ground effect are calculated with the use of 
Vortex-lattice method. Comparison between the lift-to-drag 
ratio of wings and planing body is done. The total lift force 
of the aerodynamic wings with respect to craft weight is 
determined. Finally, the idea of using the planing-wing hybrid 
craft is discussed and the conditions in which aerodynamic 


lift force of wing is beneficial for planing mono-hull crafts 
are formulated. 


Fig. 1. A schematic image of the model planing-wing hybrid craft 


THE PROBLEM STATEMENT 


This research, based on considering the issue of hybrid 
craft as a combination of aerodynamic and hydrodynamic 
forces and taking into account the Seaphantom as a model, 
is aimed at examination of usefulness of this kind of craft. 
The problem is defined as follows: a planing craft equipped 
with two wings is designed as a hybrid craft moving in water 
with the speed in the range of planing mode. The air is blown 
under the wings due to craft’s speed thus the aerodynamic lift 
is generated. A small gap between wings and water surface 
results a higher aerodynamic lift due to ground effect. Some 
further assumptions are as follows: 

- The craft has to remain in contact with water. 

- Size and geometry of the main hull follow the principles 
of planing body. 

- Dimensions and size of aerodynamic wings are in 
consistency with general arrangement of a planing craft 

- Trim angle of planing body and angle of attack of wing are 
selected rationally. 

- It should be noted that around wings water spray occurs 
that is not considered here. 

- In this study craft steady motion at constant speed in calm 
water is considered. 

- The Vortex-lattice method is used as a tool for both 
hydrodynamic and aerodynamic force calculations. 

- Friction forces are calculated by a regression formula and 
added to Vortex-lattice method calculation results. 

- Gap between water surface and the lowest point of wing 
is related to wing chord at root section, and marked H/C, 
where H is gap height and C is chord length of wing root 
section. 


THE GOVERNING EQUATION AND 
COMPUTER SOFTWARE 


Governing equation 


The Vortex-lattice methods are based on solutions of 
Laplace’s equation. They were first formulated in the late 
30s and first called “Vortex Lattice” by Faulkner in 1943. 
The velocity potential should satisfy Laplace’s equation as 
follows: 


vp = 0 (1) 


where: © is total velocity potential and, for a steady flow, 
written as follows: 


®=Ux+o (2) 


where: U is the uniform stream velocity and 9 is the perturbation 
potential. Both ® and 9 satisfy Laplace’s equation. 
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Computer software 


To perform aerodynamic calculations of a ground-effect 
craft a computer software called Autowing based on the 
boundary element method (called VLM), has been developed 
by Korneev [11]. In the computer program the main hull and 
wings are meshed by square panels on which singularities 
are distributed. This software has been developed specially 
to design wings for ground-effect crafts. The Autowing is 
capable of calculating hydrodynamic and aerodynamic lift and 
drag. The frictional drag could not be calculated by means of 
the Vortex-lattice method, hence it is computed by using an 
empirical formula inserted in the Autowing.[11]. 


THE PROCESS OF RESEARCH 


A planing body whose geometric characteristics are 
presented in Tab. 2, is selected as a hull. Its dead-rise angle 
varies from 10 degrees aft to 15 degrees amid-ship. The body 
lines of the craft and its 3D image are shown in Fig. 2 and 3, 
respectively. White area in Fig. 3 distinguishes the underwater 
part of the hull on even keel condition. 


Tab. 2. Main particulars of the hybrid craft 


Particulars 
LOA 
LBP 

B(max) 
Draught 
Height 
Displacement 
LCG from aft 


Quantity 


Fig. 2. Body lines of the hybrid craft 


Fig. 3. 3D image of the hybrid craft 
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The 3D model was used to perform hydrostatic calculations 
and consequently to obtain centre of buoyancy and 
displacement of the craft. Furthermore, its wetted surface was 
calculated at different trim angles and draughts. To determine 
hydrodynamic lift force the vessel was modelled by using the 
Autowing, as shown in Fig. 4. The calculations were done for 
four draughts and five trim angles at six different speeds as 
given in Tab. 3. 


Fig. 4. A model of the planing hull part used as input 
to the Autowing software 


Tab. 3. Input parameters of example calculations 


Quantity 
0.5, 0.65, 0.8, 0.95 
2,3, 4,5, 6 
25, 40, 55, 70, 85, 100 


Draught 


Trim 


Speed 


According to the research on the ground effect [12], 
the pressure side of wing surface, which has been carefully 
examined in the case study, should be flat. The wing shape 
is like an inverse delta, known as Lippisch type, whose 
important advantage is providing the craft with a better 
stability [13]. 

The aerodynamic calculations were also done by using the 
Autowing. The wing geometry is shown in Fig. 5. In Fig. 6 
a typical airfoil section of the wing is shown. The dimensions 
of the wings are selected proportionally to the dimensions 
of the craft. Tab. 4 shows the ratios of wing dimensions in 
respect with the craft, where C and S are the chord length and 
span of the wing, respectively. Furthermore, H is gap height, 
and L and B are the craft length and beam, respectively (see 
Fig. 7). 


Tab. 4. Wing main dimensions 


0.5, 0.75, 0.9 
0.5, 1, 1.5 


0.05, 0.059, 0.088 
8.5, 12.75, 15.3 


For aerodynamic analysis the gap of 75 cm measured from 
water level to the end point of the wing, is considered because, 
the maximum ground effect happens at a height of less than 
10 percent to the wing chord length [14]. To obtain better 
performance and achieve higher lift force the wings with the 
initial angle of attack of 7 degrees, are assumed. It means that 
if the craft is on even keel in water, its wing will have the angle 
of attack of 7 degrees. 


RESULTS OF PARAMETRIC STUDIES 


A set of calculations was done for several configurations 
of the hybrid craft, in which its draught, trim, wing size and 
speed were systematically changed. Their results are below 
presented and discussed. 


Fig. 5. Wing geometry modelled by the Autowing san a a a odd 

Since two elements support craft’s weight by dynamic lift, 
the lift-to-drag ratio is a good measure of merit for comparing 
different lifting surfaces. Obviously, the higher lift-to-drag ratio 
the better performance of a lifting surface. 

In Fig. 8 through 10, the lift-to-drag ratios are presented 
in function of angle of attack and wing area for three different 
speeds. The speeds of 25, 55 and 100 knot were selected to 
cover the ranges of low, intermediate and very high speeds. It 
is observed that the lift-to-drag ratio of wing varies from 2.5 
Fig. 7. Parameters defined for ground effect to 9.5 and decreases as the angle of attack increases. The ratio 

10.000 7 


Fig. 6. Typical airfoil section used for wing 
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Fig. 8. The lift-to-drag ratio of air wing with respect to the angle of attack and wing area at speed of 25 knot 
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Fig. 9. The lift-to-drag ratio of air wing with respect to the angle of attack and wing area at speed of 55 knot 
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Fig. 11. The lift-to-drag ratio of a planing hull with respect to speed and trim at the draught of 0. 5 m 
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Fig. 12. The lift-to-drag ratio of a planing hull with respect to speed and trim at the draught of 0. 95 m 


POLISH MARITIME RESEARCH, No 4/2012 


Ly/Ly 


La/Ly 


14.000 


12.000 


10.000 


8.000 


L/D 


6.000 


4.000 


2.000 


0.000 


20 30 40 50 70 80 90 100 


60 
V(knot) 


Fig. 13. The lift-to-drag ratio of a planning hull with respect to speed and trim at the draught of 1.36 m 
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Fig. 14. The lift ratio versus the surface ratio for different speeds at the draught of 0.5 m 
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obtains its peak value in the case of the wing of 53m? area, 
whose chord length is the least (of 50 percent of craft length) 
and span is the largest (of 150 percent of craft breadth). On 
this basis, it is concluded that the wing should be so designed 
as to obtain the possible highest aspect ratio. 

The lift-to-drag ratio of a typical planning hull is shown in 
Fig. 11 through 13, which varies from 2 to 12. It decreases as 
the draught and speed of the craft increases. Additionally, the 
greater trim the greater lift-to-drag ratio. 

A comparison of the two groups of curves leads to an 
important conclusion that there is no certain advantage from 
using the wing aerodynamic force instead of the hydrodynamic 
lifting surface force. 


Comparison of Hise ih lift with 
aerodynamic lift 


The ratio of aerodynamic lift (of wings) and hydrodynamic 
lift (of planing hull), hereafter called the lift ratio, versus the 
ratio of wing area and planing wetted surface, hereafter called 
the surface ratio, is shown in Fig. 14. 

As trim angle increases the lift ratio increases too. This is 
due to the positive angle of attack of 7 degree of the wings for 
the craft on even keel. That is why the maximum lift ratio occurs 
at the trim angle of 2 deg at the speed of 100 knot. However, 
it may not be greater than 40 percent. 

As far as the surface ratio has been concerned, increasing the 
surface ratio leads to increasing the lift ratio. Furthermore, the 
lift ratio has three distinct zones not associated with different 
trim angles. In the case of a low surface ratio, less than 1.0, the 
hydrodynamic forces plays a dominant role in supporting the 
craft weight while the aerodynamic force is negligible. For the 
surface ratio of a moderate value from the range of 1.5 up to 
2.0, the lift ratio is less than 10 percent. For the surface ratio of 
more than 2.0 the lift ratio may reach as much as 35 percent. 

(A: trim angle of 2 deg, B: trim angle of 3 deg, C: trim angle 
of 4 deg, D: trim angle of 5 deg and E: trim angle of 6 deg) 


CONCLUSIONS 


A typical hybrid craft was analyzed to examine the 
conditions in which aerodynamic lift force of wing is beneficial 
for planing mono-hull. A typical planing body fitted with wings 
of different dimensions proportional to the main dimensions of 
the craft, was considered. The analysis showed that if a planing- 
wing hybrid craft has to be used, its surface ratio should be 
more than two in conjunction with a high aspect ratio of wing. 
The main conclusions of this study are as follows: 

- Inthe planing-wing hybrid craft the hydrodynamic lift-to- 
drag ratio is usually slightly greater than the aerodynamic 
lift-to-drag ratio. Therefore, there is no comparative 
advantage from using air wing. 

- The absolute lift force developed by the hybrid craft is a key 
element for supporting its weight. Generally, the maximum 
share of the lift force produced by the air wing is equal to 
about 20 percent of the hydrodynamic lift force in practical 
conditions. 

- Some phenomena such as the hitting of spray to wings 
as well as impact of the wings against wave surface may 
produce some undesirable behaviour of a hybrid craft. 


Finally, it should be mentioned that air wing may be used 
only in special circumstances, but not in general. 

This conclusion is based on results of technical investigations 
only, but sometimes aesthetic view of a craft is more important 
than its efficiency. A planing-wing hybrid craft seems to be 
much more beautiful than a simple planing craft. 
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ABSTRACT 


In the paper [1] assumptions for fatigue life calculations in the low- and high- cycle fatigue ranges of metal 
alloys, have been formulated. Three calculation methods: that expressed in stresses, strains and a mixed 
(stress-strain) one, have been there presented. In this paper results of fatigue life calculations of C45 steel 
performed in line with the above mentioned methods have been described, and their comparative analysis 
has been made. The calculations have been carried out for two-level and multi-level load programs of the 
varying parameters: maximum load within a load program and different values of spectrum filling factor. 
The comparative analysis of the results of fatigue life calculations in compliance with the above mentioned 
methods makes it possible to assess differences in the results and depenendence of the differences on values 
of load program parameters. 


Keywords: fatigue life calculations; low- and high- cycle fatigue; C45 steel 


INTRODUCTION 


Operational loads of structural elements are generally 
random. Load spectra elaborated in accordance with relevant 
cycle-counting methods contain sets of sinusoidal cycles of 
different parameters, especially different values of the load 
amplitude S, and mean value S, [2]. In the case when the load 
Smax = Sm + S, < R, fatigue life is calculated with the use of 
the method which applies Wöhler diagram (stress approach 
— the high-cycle fatigue range, HCF), whereas for Smas > Re 
the method using Manson-Coffin diagram is applied to fatigue 
calculations (strain approach — the low-cycle fatigue range, 
LCF) [3, 4]. 

Load spectrum usually contains cycles of different share of 
Smax Values occurring both in the HCF and LCF range. Therefore 
in the authors’ research project [5] has been introduced the 
theme ofa hybrid calculation method in which fatigue damage 
due to cycles in the HCF range is calculated by means of the 
stress approach method whereas in the LCF range the strain 
approach method is used. Then the question appears whether 
applying either stress approach methods or strain approach ones 
for the entire range of loads (LCF + HCF) one obtains fatigue 
life results significantly different from those calculated by using 
the hybrid method. Significance of the so formulated problem, 
as shown in [6, 7], consists also in that criteria for qualifying the 
loads either to LCF or HCF range are ambiguous — blurred. This 
work is aimed at answering the above formulated question. 

Assumptions for calculations of fatigue life of metal alloys 
in the conditions when loading occurs both in the LCF and 
HCF ranges, were presented in the paper [1]. In the calculation 


method described in the paper three paths were distinguished: 
ls — based on application of the full (complete) Wöhler 
diagram (i.e. acc. stress approach) 2™4 — based on application 
of Manson-Coffin diagram (i.e. acc. strain approach), and 34 
— hybrid one. 

Fatigue life calculations in the varying amplitude load 
conditions require to know a load spectrum or program, 
cyclic properties of material, usually in the form of the above 
mentioned fatigue diagrams, as well as to assume an appropriate 
hypothesis for fatigue damage summation [8]. Out of many 
known hypotheses, for the calculations described in this paper 
the Palmgren-Miner linear summation hypothesis, the best 
experimentally verified one, was selected. 

In accordance with the hypothesis the sum of damages 
resulting from the number of cycles, nọ, during realization of 
loading program, is calculated from the formula: 

- for 1s path: 


yas 
b=) 1 
i=l N; -i 
- for 2"4 path: 
k 
n. 
D= 2 
i Na 
- for 3" path: ' i 
l k 
nN. Nn. 
a ae Oey 3) 
iat Ng i=(1+1) N; 


Fatigue fracture will occur when the sum of damages 
D = 1.0. which leads to the number of program repetitions: 
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d= 1.0/D, (4) 


and to the fatigue life expressed by number of cycles: 


N.=: n, (5) 
DATA FOR CALCULATIONS OF C45 STEEL 
FATIGUE LIFE 
Load programs 


The calculations were performed for load programs of two 
kinds: two-level one of variable parameters, and multi-level 
one (Fig. 1). 

In the two-level program the following parameters were 
assumed: 

- four valuess of the maximum stresses Sax Within the 
program: 650 MPa, 570 MPa, 460 MPa and 340 MPa, 

- three ratios of stress values, S,,/S,,: 0.75; 0.5; 0.25, 

- four ratios of numbers of load cycles, n01/n0: 0.75; 0.5; 

0.25; 0.1. 


In the further part of this work the assumed values of Smax 
loads are marked as follows: I — 650 MPa, II — 570 MPa, III 
— 460 MPa i IV — 340 MPa. 

Such selection of the parameters makes it possible to realize 
a wide research program comprising 64 different cases. 


a) 


b) 
Sai/ Sal 


Eac i/ Eacl 


Noi/Ny = 0,1 
€ =0,55 


Fig. 1. Schematic diagrams of the load programs assumed for calculations: 
a) two-level program, b) multi-level program 


In the tests according to 2™ path (i.e. strain approach) and 34 


path (i.e. hybrid one) values of the total strains ¢,, correspond 
to relevant values of the stresses S,. 
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The parameter which characterizes load intensity is the 


spectrum filling factor described by the formula: 
k 


C= >a Bat (6) 
i-l S a max n o 

which, on transformation, for the two-level program, takes the 

following form: 


S 
c= ta „Sifi Bal (7) 


ng Sai n 


Values of the spectrum filling factor for the assumed 
parameters are collected in Tab. 1. 


Tab. 1. Values of the spectrum filling factor € 


The data contained in Tab. 1 indicate the wide range of 
loading conditions for fatigue calculations, starting from 
¢ = 0.325 which corresponds to e.g. the loading of track system 
elements of road vehicles, up to ¢ = 1.0 which corresponds to 
sinusoidal load of constant amplitude. 

The multi-level load program was elaborated on the basis 
of the load spectrum of the mean value of the spectrum filling 
factor C, calculated from the corresponding data given in Tab. 
1. And, the graphical form of the program is shown in Fig. 2 
in the system of relative quantities. 

The program assumed for the tests in question has 10 load 
levels (S,; or £i) of the same cycle capacity n,/n, = 0.1 and 
the same load span S,-S,;,,= 0.1 © S,,. The spectrum filling 
factor 6 = 0.55. 

Moreover, for the fatigue life tests and calculations were 
assumed maximum stress levels complying with those taken for 
fatigue life analyses in the two-level loading conditions. 


ait] 


Static and cyclic properties of C45 steel 


Mechanical properties of C45 steel were determined on 
the basis of specimens prepared in compliance with the Polish 
standards: static properties — acc. PN-EN ISO 6892-1:2010. 
cyclic properties — acc. PN-84/H-04334. The values of the 
determined parameters are: R,, = 682 MPa, R, = 458 MPa, E= 
2.15:105 MPa. The Wöhler diagram experimentally determined 
for R = -1.0 is described by the formula: 


logS, =-0.1020logN + 2.9611 (8) 


for which the fatigue limit at N, = 106 is equal to Sg = 223.5 
MPa, and the exponent m, = 9.8. 

The Manson-Coffin diagram experimentally determined is 
described by the formula: 


1204 

= 5 (2 N r) 
2.15-10 

To transform the load program from that formulated in stress 


approach to strain approach one, was used the experimentally 
determined Ramberg-Osgood diagram of the following form: 


5.06 
_ Ss (Sa (10) 
Eac = 3 + 
215-10° (1232 


—0.1033 — 0.475 


+0.2179QN,) 4” 0) 


ac 


RESULTS OF CALCULATIONS 
FOR TWO-LEVEL LOADS 


Results of calculations according the 1%“ path 


The calculation results in the form of Gassner diagrams are 
shown in Fig. 2 on the Wohler diagram background. As results 
from the data assumed for the calculations, the levels I, II and 
II correspond to the LCF condition, i.e. S,, > R,, similarly the 
values S,, for S,,/S,, = 0.75 are contained within the low-cycle 
fatigue (LCF). The remaining values, S,, and S,, are smaller 
than R,, hence in compliance with the assumed criterion they 
are rated to be in the high-cycle fatigue range (HCF). 


Stress amplitude S,, MPa 
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----- = --- - - - - - +--+ 
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Fig. 2. Fatigue life diagrams for C45 steel determined by calculations 
according to 1s path in the two-level loading conditions: 1) Wöhler diagram 
for n/n = 1.0; 2) Gassner diagram for n/n = 0.75; 3) Gassner diagram 
for n/n = 0.5; 4) Gassner diagram for n/n, = 0.25; 5) Gassner diagram 

for n/n = 0.1 


Results of calculations according to the 2"4 path 


In compliance with the description given in p. 1. in the 2.4 
path calculations the load ranges LCF and HCF have to be 
related to Masnon-Coffin diagram. To this end, it is necessary, 
for description cyclic strain diagram, to transform loads 
expressed in stress units into those expressed in total strains by 
using the Ramberg-Osgood diagram (10). The total strain values 
corresponding to relevant stress values are given on the ordinate 
axis of the diagrams in Fig. 3. The number of 2N; recurrence 
to fatigue fracture was calculated by using the Mason-Coffin 
formula (9). The calculation results are graphically presented 
in Fig. 3. 

In Fig. 2 and 3 the diagrams are presented in the form 
of shadowed bands. Their left edges correspond to the ratio 
S/S, = 0.75, right ones — the ratio S,,/S,, = 0.25, and the 
middle line of the band corresponds to the ratio S/S, = 0.5. 


Calculations in compliance with the hybrid 
method for two-level loading (the 3” path) 


According to the hybrid method (the 3% path) complete 
damage corresponds to sum of the damage resulting from the 
loading within the low-cycle range (LCF), calculated by using 
Manson-Coffin fatigue diagram and the damage due to the 
loading within the high-cycle range (HCF) calculated by means 
of Wöhler fatigue diagram according to the formula (3). 

The sums of damages calculated acc. (3) and the fatigue 
life calculated in compliance with the hybrid method, are 
presented in Tab. 2. 


Strain £ac 
10° Mterte meet eee mme mremmeeemmemmmmmmmermmmmemmeemeeemem peenema prmmemmemeemmmeeempeemeemme meenemen 
3 i i 5 i p 
Level I (Eac1 = 4,24-10°) 
10° 
10! 10° 10° 10° 10° 10° 10’ 
Number of cycles N 


Fig. 3. Fatigue life diagrams for C45 steel determined by calculations 
according to 2"4 path in the two-level loading conditions: I) Manson-Coffin 
diagram; 2) Gassner diagrams for n/n, = 0.75; 3) Gassner diagrams for 
n/N, = 0.5; 4) Gassner diagrams for n/n, = 0.25; 5) Gassner diagrams 
for n/n = 0.1 


In line with expectations, for high load values (the level I) 
the results of fatigue life calculated by using the hybrid method 
(the 3 path) are closer to the results of calculations according 
to the 2™4 path, and for low load values (the level IV) the results 
of fatigue life calculations according to the 34 path are closer 
to those calculated in compliance with the 1st path. 

The above given statement shows that in the case of higher 
values of the ratio the higher level, i.e. the amplitude S,, or €,,, is 
significantly decisive of fatigue whereas influence of the lower 
level, i.e. S,, or €,,.. is not significant. This is confirmed by the 
data given in Tab. 3 and 4. Tab. 3 contains the numerical data 
concerning the ratio of the fatigue life calculated by applying 
stress approach (the 1s‘ path) and that calculated by using strain 
approach (the 24 path). 

As results from the data contained in Tab. 3, for high stress 
values and strains corresponding to them (the load levels I and 
II) the more conservative results are achieved from the strain 
approach calculations (the 274 path), whereas for lower load 
values (the levels II and IV) the results of stress approach 
calculations (the 1s path) are more conservative. Fig. 4 
illustrates the conclusion. 


2.0; N-SN-E 


i 
è - 0.5; 


0 01 02 03 04 05 06 07 08 09 10 


Symbol S/S e - 0.75; @ - 0.25 


Fig. 4. Dependence of the ratio of the fatigue life calculated acc. the 1" path 
and that calculated acc. the 24 path on load level and the test program form 
described by the ratios S,,/S,, and ng /ng 


As results from the diagrams shown in Fig. 4, the influence 
of loading program form on difference in the calculation 
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level 


Tab. 2. Values of damage parameter and fatigue life calculated in compliance with the hybrid method (the 3 path) 


sem 


0.433540 
0.428100 
0.428000 


0.084850 
0.083423 
0.083400 


0.297120 
0.286198 
0.286002 


0.058500 
0.055655 
0.055600 


0.144466 
0.143098 
0.142800 


0.032260 
0.027882 
0.027800 


0.085180 
0.083434 
0.083400 


0.059100 
0.055667 
0.055600 


1692.0 
1796.4 
1800.0 


0.032950 
0.027701 
0.027600 


0.077200 
0.057557 
0.572004 


0.016330 
0.011198 
0.011100 


0.011743 
0.011112 
0.011100 


8515.7 
8999.1 
9009.0 


4.76:10° 
4.69-10° 
4.69:10° 


30166 
31912 
31912 
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1.667-10% 
1.667-10% 


6.281-10% 
6.250-10% 
6.250-10% 


159200 
160000 
160000 
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Tab. 4. The ratio of the fatigue life calculated acc. stress approach (the 1s path) and acc. strain approach (the 2"4 path), 
respectively, and the fatigue life calculated by using the hybrid method 


Ny,/Ny 


S.2/Sa 


No. of level 
Ead Exi 


1.0 


0.75 


0.5 0.25 0.1 


NAN. 


NZ/NU 


NS/NU 


No/NEINZ/Nid|NS/NEINZ/Nd|NS/NUINZ/N¢ 


4.24 - 102 


2.22 + 102 


8.98 - 10-3 


3.06 - 103 


results according to the 1st path and 24 path, expressed by the 
fatigue life ratio (NŠ/N®), is not significant and observed only 
for the programs of ny,/ny = 0.1. However in accordance with 
expectations, the load level influence on the ratio (N3/N®) 
is significant and amounts to: (NS/N®) z 1.6 for the level I, 
(NS/N®) = 1.4 for the level II, (N3/N®) = 0.92 for the level III 
and (N/N?) = 0.71 for the level IV. 

On assumption that the reference point is the fatigue life 
value NË, calculated in accordance with the hybrid method 
(the 3* path), the fatigue life ratios NS/N'! and NẸ/N} were 
calculated. The data are collected in Tab. 4 and illustrated by 
the diagram shown in Fig. 5a and 5b. 

Like in the case of the data shown in Fig. 4, load level shows 
significant influence on the conformity between the results of 
calculations according to the 1%, 274 and 3" path, and a lower 
influence is produced by load program form described by the 
ratios S,,/S,, and no,/no_ the latter - with the exception for no,/ 
ny = 0.1. Comapring the diagrams in Fig. 5a with the relevant 
diagrams in Fig. 5b one states that a greater conformity is shown 
by the calculation results according to the 1st and 3 path for 
low load levels (obviously, full conformity is observed for the 
load levels III and IV) as well as the 2»4 and 3" path for high 
load levels (obviously, full conformity is observed for the load 
levels I, II and HI). 

The maximum difference between fatigue life calculated 
according to the 1stand 34 path amounts to 6 = 39%, whereas 
the relative difference between fatigue life calculated according 
to the 2™ and 3 path amounts to 6 = 33%. 


CALCULATIONS FOR C45 STEEL IN 
MULTI-LEVEL LOAD CONDITIONS 


Results of calculations according to the I* path 


The multi-level load program has been already described 
in p. 2.1. In Fig. 6 is presented Gassner diagram (2) on the 
background of Wöhler diagram (1), moreover for comparison, 
by means of the dotted line (3) is shown the Gassner diagram 


a l N 
Symbol S2/Su: @ -0.75; @-0.5; @-0.25 


Fig. 5. Dependence of the ratio of the fatigue life: a) calculated acc. the 1st 
path and that calculated acc. the hybrid method (the 3% path), b) calculated 
acc. the 24 path and that calculated acc. the hybrid method (the 34 path), 
on load level (I, II, III and IV, respectively) and the test program form 
described by the ratios S/S, and n/n 


for the same two-level loads and the same spectrum filling 
factor C,. 
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As results from the comparison the fatigue life determined 
in the two-level load conditions differs from that determined 
in the multi-level load conditions within the entire variability 
range of S,- 


Stress amplitude Sa, MPa 


10! 10° 10° 108 10° 10° 107 
Number of cycles N 


Fig. 6. Fatigue diagrams for C45 steel in stress approach: 1) Wöhler 
diagram, 2) fatigue life diagram (acc.Gassner) for the multi-level load 
program of the spectrum filling factor ¢ = 0.55, 3) fatigue life diagram 

(acc.Gassner) for the two-level load program 
of the spectrum filling factor € = 0.55 


In accordance with the criterion of LCF assumed in the 
program for the load level I (S,, = 650 MPa) four first levels 
of the program are situated within LCF range. Similarly, in the 
program for the load level II (S,, = 570 MPa) three first levels 
of the program are situated within LCF range, whereas in the 
program for the load level III (S,, = 460 MPa) only the first 
level is contained within LCF range. The remaining levels in 
all the programs are situated in HCF range. 


Results of calculations according to the 2"¢path 


The Gassner fatigue life diagram (the line marked 2) on the 
background of the Manson-Coffin fatigue diagram (the line 
marked 1) is presented in Fig. 7. On the discussed diagram 
is presented for comparison the fatigue life diagram for two- 
level load program (the line marked 3). Both the two-level 
and multi-level program are of the same value of the spectrum 
filling factor ¢. The mutual location of the diagrams 2 and 3 
shows small differences in fatigue life values calculated for 
multi-level and two-level program. 

Like in the case of calculations according to the 1s path, 
some part of the programs is located in LCF range and the 
remaining in HCF range. Specification of the levels has been 
given in the discription of the 1st path. 


Results of calculations according to the 3" path 


The hybrid method (the 3" path) consists in summing 
damages in LCF range in accordance with the method based on 
the concept of application of Manson-Coffin fatigue diagram 
(the 24 path), and in HCF range in accordance of the method 
based on the concept of application of Wöhler diagram (the 
1stpath). 

In Fig. 8 is presented the comparison of fatigue life diagrams 
for C45 steel, determined in the conditions of multi-level load 
program. The diagrams were elaborated in the system: the 
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Fig. 7. The fatigue life diagram for C45 steel in the conditions of the 

multi-level load program (line 2) on the background of the Manson-Coffin 


fatigue diagram (line 1), as well as the diagram for the two-level load 
program (line 3) 


load level S, (1, IL, Ill, IV) versus fatigue life expressed by the 
number of cycles N,. As results from the comparison of the 
diagrams, the results of the stress approach calculations (the 
1s path) are closest to the values of fatigue life determined by 
using the hybrid method (the 3‘path). Only for the level I (650 
MPa) the results of the strain approach calculations (the 24 
path) show a higher conformity with those calculated by using 
the hybrid method (the 31 path). The above given observation 
is in compliance with the conclusion formulated in the case 
of analysis of the data achieved from the calculations in the 
conditions of two-level load programs with the exception that in 
the case of multi-level load programs differences in calculation 
results are much smaller and amount, in extreme cases, to: 

- forthe level I (S,, = 650 MPa), between the results acc. the 

34and 1st path: 


- forthe level I (S,, = 650 MPa), between the results acc. the 
34¢and 24 path: 


- for the level IV (S,, = 340 MPa), between the results acc. 
the 34and 1s path: 


- for the level IV (S,, = 340 MPa), between the results acc. 
the 34and 24 path: 


ANALYSIS OF RESULTS AND SUMMARY 


The basic analysis of fatigue life calculation methods for C45 
steel concerns programmed two-level loads whose parameters 
assumed in the program: S,,; S,./S,, and no,/ny made it possible 
to perform calculations for 64 cases. The so wide range of the 
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Fig. 8. Gassner fatigue life diagrams for C45 steel in the multi-level load 

conditions, determined with the use of the following calculation methods: 


1) acc. stress approach (the 1* path) 2) acc. strain approach (the 2"4 path), 
3) hybrid one (the 3” path) 


two-level load cases makes it possible to quantitatively assess 
influence of program’s parameters on fatigue life calculated in 
accordance with all the calculation paths. 
As expected, results of the calculations demonstrate: 
- prevailing influence of the load level: S,, or €,,; 
- much lower influence of the program’s form: S,,/S,, and 
Ny,/Ny Which decreases along with spectrum filling factor 
increasing. 


Purposefulness of assuming, in the calculation program, the 
loads of high values of the spectrum filling factor ¢ (Tab. 1) 
results from the necessity of stressing significant contribution 
of loads from LCF range. This is the case which corresponds 
to trends in the proposed calculation methods. The programs 
in question differ, as to the range of C -values, from the known 
programs elaborated on the basis of random operational loads 
of structural elements. For instance [9] for the track system 
elements of road vehicles the ¢ factor takes values from the 
interval of 0.25 + 0.35. 

As results from the data contained in Fig. 2 (the 1* path), 
Fig. 3 (the 2™4path) and Tab. 2 (the 3" path), influence of fatigue 
damages due to cycles of the 2»4 degree amplitude (S,,; €,..) on 
fatigue life is insignificant, with the exception of the program 
of no)/N) = 0.1. 

As expected, is observed a high conformity of the calculation 
results according to the assumed paths on the load level 
Sai = R, (Fig. 4 and 5) for which the fatigue life ratios NS/N€, 
NS/N# and N&/Ni are close to 1. 

The wide range of analysis of calculation results of fatigue 
life for the two-level load programs made it possible to 
elaborate an appropriate multi-level load program (of 10 levels). 
The selected ten-level load program (Fig. 1b) is characteristic 
of the same number of cycles on particular levels (no; = n,/10) 
and the same load increase on particular levels (AS, = S,- 
Saan = Sa1/10). It leads to the spectrum filling factor € = 0.55, 
close, as for value, to the mean counted from Tab. 1 for two- 
level load programs. Moreover the same load levels S,, and ¢,,,. 
were maintained, that makes it possible to conduct comparative 
analysis of calculation results of fatigue life for multi-level 
programs and those for two-level programs. 

As results from the data contained in Fig. 6, the low values 
of the stress amplitude S,, corresponding to the program levels 
from i = 6 through i = 10, are of insignificant influence on the 


summing of fatigue damages (less than 0.001 D,,) and therefore 
may be neglected in fatigue life analyses. The above formulated 
statement is highly important for programmed fatigue tests 
since exclusion, from the program, of load levels which 
insignificantly affect fatigue life, shortens testing duration 
time considerably. 

The similar conclusion results from the calculations 
according to the 24 path (Fig. 7) and those according to the 
3 path (Fig.8). From comparison of the fatigue life diagrams 
determined in accordance with the analyzed paths it results, 
like in the case of two-level load programs, that the results 
of the strain approach calculations (the 24 path) for high 
load levels (the level I) are closest to those obtained from the 
hybrid method (the 3" path) whereas for low load levels (the 
level IV) the results of the stress approach calculations (the 1st 
path) are closest to the results of the calculations according to 
the hybrid method (the 3rd path). For mean values of loads 
(the level II and HI) the results of the fatgue life calculations 
according to the assumed paths (1%, 2»4 and 34 one) are close 
to each other. The relative maximum differences in calculation 
results amount, in extreme cases, to: between those acc. the 34 
and 1“ path — the relative difference 81, ,=—0.42; between those 
acc. the 3 and 2» path — the relative difference 51), = —0.66, 
and are significantly smaller than the corresponding differences 
observed in two-level programs. 

In Fig. 6 are collected fatigue life diagrams determined 
by means of calculations in accordance with the 1* path 
for the multi-level load program (the line 2) and two-level 
program (the line 3) of similar values of the spectrum filling 
factor ¢. From their comparison rather low relative difference 
in calculated fatigue life values, amounting to 0.40 in the 
extreme case (for the load level IV), results. In the case of the 
calculations according to the 24 path relative difference in 
calculated fatigue life values (for the load level IV), amounts 
to 0.26. The above formulated statements indicate that to 
simplify load program form by reducing number of program 
levels is possible, that is very important for programmed 
experimental tests. 

The problems described in this work find application 
to fatigue calculations of ship structures and significantly 
widen the calculation methods described in the subject-matter 
literature [10, 11 and 12] to cover the load range common for 
both low-cycle and high-cycle fatigue. The above mentioned 
load case - a sum of regularly changeable loads and high 
overloads imposed on them - is that commonly occurring in 
ship structures. 
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NOMENCLATURE OF MAJOR NOTATIONS 


— fatigue damage due to realization of n, cycles of loading, 
— number of cycles to fatigue damage — general notation 
(fatigue life), 

— fatigue life expressed by number of cycles, determined in 
conditions of programmed loading, 

number of cycles to fatigue fracture read from Manson- 
Coffin fatigue diagram for the total strain £i» 

— number of cycles to fatigue fracture read from Wöhler 
diagram for the stress amplitude S; 

— basic number of cycles corresponding to fatigue limit, 

— general notation of stress, [MPa], 

maximum stress value in sinusoidal cycle, [MPa], 

— minimum stress value in sinusoidal cycle, [MPa], 
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S — stress amplitude in sinusoidal cycle, (S, = 0.5(S nax Smin)» 
[MPa], 

S — mean stress value in sinusoidal cycle, (S, = 0.5(Sija-Sinin))» 
[MPal, 

R — stress ratio (R = S nin Smax), 

R,  — yield point [MPa], 

R — tensile strength of a material [MPa], 

k — number of load levels in a loading program, 

l — number of load levels in a loading program of LCF range, 


n  — number of cycles in a loading program, 
No,  — number of cycles on 1* level of a loading program, 
Ny;  — number of cycles on i-th level of a loading program, 
Ma — exponent in Wöhler diagram formula, 
— relative difference between results of tests and 
calculations, 
£ — total strain value, 
C — spectrum filling factor, 
A — number of repetitions of a program to fatigue fracture, 


LCF — low-cycle fatigue, 
HCF — high-cycle fatigue. 
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ABSTRACT 


The fluid-solid interaction (FSI) dynamic responses for a Light Weight Pyramidal Sandwich Plate Structure 
(LWPSPS) under different water-entry velocities (Im/s-6m/s) are investigated numerically and theoretically. 
The characteristics of impact pressure and structure deformation are obtained by using LS-DYNA code 
based on the proposed 3D multi-physics (air-water-solid) FEM model. Numerical results show that the 
average water impact pressure of LWPSPS is much lower than that of the monolithic plate with same mass. 
Moreover, a phenomenon called “local air cushion” is observed for LWPSPS which does not exist for 
a monolithic plate. Theoretical hydroelasticity analysis demonstrates that the FSI effect of LWPSPS is weak 
when the ratio of water impact duration time to wet natural vibration period is greater than 4. In the study, 
an engineering estimation method to predict the maximum deformation of the LWPSPS is proposed, in which 
the total deformation is divided into two parts, i.e. local field deformation and global field deformation, 
and they are both computed using analytical model. Good agreement between the numerical results and 
ones obtained from the proposed engineering estimation method is achieved. Furthermore, the geometric 
variation sensitiveness analysis is also conducted. 


Keywords: light weight pyramidal sandwich plate structures (LWPSPS); multi-physics numerical model; 
maximum deformation; quasi-static; local air cushion 


INTRODUCTION 


The water impact phenomenon exists in many practical 
engineering areas such as slamming of ship experienced 
boisterous weather and heavy sea, seaplane landing on water 
in an emergency situation et al. [1]. The significant local 
structural deformation may occur due to high peak water 
impact transient shock pressure in a very short duration [2]. 
Because of the severity and significance of slamming accident 
on naval vessels, many research efforts have been devoted 
during the past century. Since the late-1920s, the pioneer works 
of hull-water impacts have been investigated by analyzing 
a two-dimensional wedge impact on a calm water surface. Two 
of these early studies related to slamming problem of a rigid 
v-shaped wedge with small deadrise angle were carried out 
by Von Karman [3] and Wagner [4]. After then Chuang [5] 
performed the systematic tests with flat bottom panels in free 
fall against the free surface of the water in order to overcome 
the problem of the Von Karman asymptotic theory due to the 
lack of application for zero deadrise angles and to detect the 
effect of trapped air between the free surface and the plate on 
the dynamic response of the plate. He concluded that the effects 


of air cushion cannot be neglected if the deadrise angle was 
less than 3 degree in his experiments. 

The equivalent rigid/quasi-static method has been widely 
adopted in the design rule for water impact problem due to its 
convenience [2]. However, the hydroelastic interaction effects 
should not be neglected in the design of high-speed craft [6]. 
Thus the hydroelastic criterion and models have been developed 
from 90s in [7-10]. Generally, the hydroelastic effects of 
water impact event are typically characterized by a relation 
between the loading period and the natural vibration period of 
structure [7-10]. Most of the previous hydroelastic interaction 
investigations are limited to simple steel structures such as 2D 
transverse section with v-shaped section of a ship-hull. 

In recent years, the mechanical behavior of composite 
materials and sandwich structure are widely concerned [11]. 
Kozak [12] presented a deep and extensive view on selected 
problems concerning application of steel sandwich plates 
in shipbuilding industry. To investigate the fluid-structure 
interaction (FSI) during the water impact process, Qin and Batra 
[13] developed a hydroelastic model based on a (3, 2)-order 
sandwich theory and 2D Wagner’s water impact theory. The 
hydroelastic effect on hull’s deflections and the unsteady water 
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Fig. 1. Light weight sandwich plate structures with pyramidal lattice cores (LWPSPS) 


impact load has been delineated. Baral et al. [14] presented 
results from atest simulating the water impact (slamming) 
loading of sandwich boat structures. Test results suggested that 
a novel foam core reinforced in the thickness direction with 
pultruded carbon fiber pins would offer significant resistance 
to water impact load. 

The lattice structures with various core topologies have now 
been recognized to be attractive candidates for multifunctional 
ultra-light structures due to their superior properties of light 
weight, high specific stiffness and excellent impulsive resistant 
performance et al. [15-18]. As the lattice structures are periodic 
micro-architectured cellular solids in which the individual 
truss elements deform by stretching rather than by bending. 
Consequently, the dynamic behavior of sandwich plates with 
pyramidal lattice truss topology cores (LWPSPS) is quite 
different from other traditional constructions, including those 
with honeycomb core or corrugated (prismatic) core et al. Most 
of the previous studies focused on the static mechanism and 
energy absorption behavior of LWPSPS subjected to intense 
impact loading [15-18]. But the fluid-solid interaction (FSD 
characteristics of LWPSPS subjected to water impact load are 
neglected, which should also be considered due to the fact that 
LWPSPS used in the construction of naval structures would 
be impacted in slamming conditions. Therefore, investigation 
should be undertaken on analyzing the hydroelastic pressure 
and deformation of LWPSPS subjected to slamming load at 
present. 

To this end, the present study first primarily concerns on 
the detailed dynamic characteristic and hydroelastic mechanism 
of sandwich plates with pyramidal lattice truss topology cores 
(LWPSPS) subjected to water impact loading. The nonlinear 
explicit finite element (FE) software LS-DYNA is employed to 
build up the 3D multi-physics (air-water-solid) numerical model 
with the water modeled as an inviscid, compressible and non- 
irrotational fluid. In order to verify the FEM model, the central 
impact pressure values of circular solid plate are compared 
with experimental results. It is found that hydroelastic water 
impact pressure predicted by LS-DYNA agrees well with those 
experimental data. Furthermore, the detailed characteristics 
of hydroelastic water impact pressure and hull deformation 
of LWPSPS are summarized and compared with those of the 
monolithic plates with equal mass. Then the hydroelastic effect 
of LWPSPS is discussed theoretically by using the criteria of 
2t/T (the ratio of water impact duration time 21, to the wet 
natural period T). Subsequent to calculate the maximum total 
elastic deformation of LWPSPS, an engineering theoretical 


32 POLISH MARITIME RESEARCH, No 4/2012 


estimation method is proposed, in which the total deformation 
is divided into two parts, i.e. local field deformation and global 
field deformation. Results from this proposed model agree well 
with those from the 3D FSI finite element analysis. Moreover, 
the individual effects of design parameters of LWPSPS such 
as face sheet thickness, thickness of lattice core and relative 
density of core on the water impact pressure and structural 
deformation are also investigated in the study. 


DESCRIPTION OF THE LWPSPS 
GEOMETRY MODEL 


To analyze the characteristic of LWPSPS, the first step is to 
choose its different constituents. LWPSPS consists of two thin 
face sheets attached to both sides of a light weight lattice core 
with an exposed area 595.1 mm x 270.5 mm referred a practical 
monolithic flat panel of cruiser [19] (see Fig. 1). The mass of 
LWPSPS (case 1 - case 6) is equal to the monolithic flat panel 
which the thickness of this mild steel monolithic flat panel is 
equal to 11.11 mm.The boundary condition of the panel can be 
taken as clamped [19]. Therefore, a sealed off round structure 
made by rigid materials (rigid out wall, see Fig. 1) is used 
to model clamped boundary condition (also see Fig. 1). The 
calculated parameters for all 22 cases are shown in Table 1. 


NUMERICAL SIMULATION MODEL 
OF LWPSPS 


The hydroelastic performance of LWPSPS subjected to 
water impact load is investigated using 3D multi-purpose, 
explicit and implicit finite element program LS-DYNA. Both 
Lagrangian and Eulerian solvers are available in LS-DYNA 
code to enable modeling of structures and fluids in a single 
model and to simulate the interaction between them. Interaction 
between the fluids and structures is achieved through arbitrary 
Lagrangian-Eulerian (ALE) description. In the fluid-structure 
coupling algorithm (typically, penalty coupling algorithm), 
two superimposed meshes are considered, a fixed Eulerian 
or ALE mesh for the fluid and a deformable Lagrangian 
mesh for the structure. Additionally, the ALE Multi-Material 
formulation is a method allowing the finite element mesh to 
move independently from the material flow and where each 
element in the mesh can contain a mixture of two or more 
different materials such as air and water. Therefore, ALE Multi- 
material formulation and penalty coupling algorithm [20-22] 
are adopted in the present 3D multi-physics numerical model 


Tab. 1. Various computational cases considered 


a(cm) 


t.xt(mmxmm) 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 
6.00x6.00 


6.00x6.00 


6.00x6.00 
6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


6.00x6.00 


5.00x5.00 


5.50x5.50 


6.50x6.50 


to simulate water impact phenomenon. The detailed description 
of ALE formulation and penalty coupling algorithm is provided 
in References [20-21]. 


Constitutions of solid 


The LWPSPS is made of mild steel for marine industry [23]. 
For the constitutive law of this steel, the kinematic hardening 
tule and elastic-linear work-hardening model are adopted in 
the analysis. And the materials properties are described as 
follows [19]: 

- Modulus of elasticity E = 206.85 GPa; 
- Yielding strength oy = 206.85 MPa; 

- Density p, = 7850 kg/m3; 

- Poisson ratio v = 0.3; 

- Tangential modulus E,= 250 MPa. 


Generally, the effect of strain rate is sufficiently important 
to analyze the dynamic responses of structures. The plastic 
kinematic hardening description is commonly used in these 
impact analyses, which is a strain rate dependent elastic—plastic 


7.00*7.00 


model. In the analysis, strain rate is taken into account by using 
the Cowper-Symonds model [23] which scales the yield stress 
by the strain rate dependent factor. According to Ref. [24], 
the strain-rate sensitivity model provides a reasonably good 
representation of test results up to strain rate of 1000 S-1. Eq. (1) 
shows the relationship between the stress ratio and the plastic 
strain rate [23]. The dynamic yield stress cyp is a function of 
plastic strain rate and can be expressed as follow: 


i ‘ 1/q 
ee ofin (1) 


Oy cs 


where the Cowper—Symonds strain rate parameters D, and q 
have values of 40.4 S- and 5 for mild steel in this study [23], 
and oy is the initial static yield stress of the material. 


Equation of state 
In general, the equation of state (EOS) specified for water is 


a Gruneisen equation and it can be used to calculate the internal 
characteristics of seawater. It properly handles wave propagation 
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Fig. 2. 3D multi-physics finite element model of air, water and LWPSPS 


phenomenon in water impact events by incorporating a non- 
linear shock velocity-particle velocity relationship. And this 
equation addresses fluid under tension/compression which can 
capture cavitation phenomena in seawater. The Gruneisen EOS 
is defined as follow: 


past t+ -t )u-(3 
Se ey 


2 
Pater ~ i 2 7 3 


=f" 
(+u) (+u) 


1-(S,-Dp-S, 


+(y+nwe 


where e is the internal energy per unit volume, C is the intercept 
of the us-up curve, S,, S,, and S, are the coefficients of the slope 
of the us-up curve, y, is the Gruneisen gamma, and 7 is the 
first-order volume correction to yọ. The constants C, S}, S}, S;, 
Yo, and y et al. are all input constant parameters [20]: reference 
sea water density p,,. = 1025 kg/m3; sea water temperature 
T = 20°C; speed of sound in water C, = 1480 m/s; dynamic 
viscosity coefficient v, = 1.13E-03 N-s/m2; intercept of us—up 
curve C = 2417 m/s; Gruneisen coefficient S, = 1.41, S, = 0, 
S, =0; Gruneisen gamma y, = 1.0; first-order volume correction 
to gamma n = 0; internal energy of water per specific volume 
e = 1.89 E + 06 J/m3; related volume at time 0 to reference 
specific volume is equal to 1. 


And the compression of sea water is defined in terms of 
relative volume V as: 


oar (3) 
Pwo 

where p„is the reference sea water density; p, is the overall 

sea water density. 


In LS-DYNA, the perfect gas law (gamma law) is used 
as the equation of state for air. This equation of state is given 
by [20]: 


P air0 (4) 

where p,;, is the overall air density for perfect gas, e, is the 
internal energy of air per unit mass and y is the adiabatic 
expansion coefficient. Here, y = 1.40, e, = 2.11 E + 06 J/kg 
and P iro = 1.239 kg/m? are conducted in the present study for 
perfect gas. 


Finite element model and verification 


Fig. 2 shows the 3D multi-physics finite element model 
adopted in the simulation. The face sheets of LWPSPS are 
modeled as Lagrangian Belytschko-Tsay shell elements (see 
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Fig. 2), where the lightweight lattice cores are modeled as 
Hughes-Liu beam elements, seawater and air are modeled as 
Euler solid elements (see Fig. 2). The dimensions of the water 
region for the 3D finite element model are 2.6 x 1.3 x 0.8 m3 
while those of the air region are also 2.6 x 1.3 x 0.8 m3. 

In general, it is known that the results based on ALE 
algorithm are sensitive to the Euler mesh density. The Euler 
mesh needs to be fine enough to capture the highest gradients in 
the pressure fields, yet a coarser mesh is favorable in terms of 
computational cost. Furthermore, for the selection of the contact 
stiffness in penalty based contact algorithm, it is required that 
the maximum pressures are approximately known ahead. So the 
non-physical contact penetration can be controlled [20-22]. 

Hence, the Euler mesh density and parameters of contact 
penetration stiffness in simulation are firstly verified by water 
impact experimental results provided in literature [25]. With 
four different impact velocities, a circular monolithic plate 
(made by 6061-T6 aluminum alloy, thickness = 11.5 mm, 
diameter = 200.0 mm) free falls on to the water surface. The 
simulation model of this circular plate is similar to that of 
LWPSPS, as shown in Fig. 2. The other detailed description of 
these water impact experiments can be found in reference [25]. 
The material constants for 6061-T6 aluminum alloy conducted 
in the numerical validation analysis of this circular plate model 
are: E = 68.94 GPa, E, = 0.7 GPa, v = 0.33, p,, = 2700 kg/m3 
and oy = 352 MPa [26]. 

Water impact velocities employed in these experiments 
are 1.40 m/s, 1.98 m/s, 2.43 m/s and 2.80 m/s respectively. 
When the Euler mesh size is about 1.0 cm, the high frequent 
oscillations in the pressure curve are not obvious, and the 
peak FSI pressure results are found acceptable compared with 
experimental results in Table 2. 


Tab. 2. Tabulated central hydroelastic pressure results 
for aluminum alloy plates 


Average 


Pressure- 
Experiment 
[KPa][25] 


Pressure- 
Numerical 
[KPa] 


62.07 | 76.66 


Relative Error | 9.08% | 0.18% 


Water impact FSI pressure 


The empirical expressions are commonly used to describe 
the hydrodynamic impact pressure. Yuzuru (1954) [25], Chuang 


(1970) [5] and Jones (1973) [27] all proposed typical equations 
which are as follows: 
1. Yuzuru (1954) [25] 


i. = KV (5a) 
2. Chuang (1970) [5] 
Pmax =A = By (5b) 
3. Jones (1973) [27] 
Pmax = kv? (5c) 


Notations A, B and k in Eqs. (5) are constants based on 
the geometry dimensions and material properties of structure. 
It can be observed that the main difference of Eqs. (5) is the 
selection of the powers of velocity v. The unified expression 
of Eqs. (5) is below: 


Pmax ~ Po = kve (6a) 


or: 


MO max — Po) = In(k) + aln(v) (6b) 


where p,, is atmosphere pressure; a and k are also the parameters 
based the geometry dimensions and material properties of 
structure. 
The water impact loading pressure of LWPSPS in terms of 
Ln (v) for six different free fall velocities (C1-C6) are shown in 
Fig. 3. For comparison, the numerical results of monolithic flat 
panels with the same mass are also given in Fig. 3. According 
to the simulation results plotted in Fig. 3, the following 
observations can be drawn: 
- The performance of central water impact pressure for both 
monolithic plate and LWPSPS follows a linear trend, as 
determined from a least square linear fit: 


In(Pax — Po) = 3.75 + 1.641n(v) 
In(Pinax — Pa) = 4.05 + 1.44In(v) (7b) 


- By means of the least square linear fit, the average 
FSI pressure for monolithic plate and LWPSPS is also 
determined: 


Inma — Po) = 2-71 +2.03In(v) (19) 


nPar — P..) = 2.63 + 1.69In(v) (74) 


- During the hull-water impact, the advantage of LWPSPS 
deceasing peak pressure amplitude on entire FSI surface 
at high impact velocity (>2m/s) is observed. But the 
difference of water impact pressure characteristics between 
LWPSPS and monolithic plate located in the central point 
is insignificant. 

- The values of a and k are constants though the range of 
water impact pressure varies widely. This conclusion has 
been identified by previous experimental results of solid 
plate [5, 25, 27]. Furthermore, the value of a of average 
FSI pressures for LWPSPS is close to 1.69, which is lower 
than the one of monolithic plate ( = 2.03). 


(7a) 


As asummary, the impulsive shapes are shown in Fig. 4. 
A more precise two parameters parabolic mathematical 
expression is defined as follows, which is different from the 
expressions conducted in the literature [27-29]: 


K,t?+K,t, O<t<2r, 


(= 
P 0, w, St 


(8a) 


where: 


p 
K=- "a K, = 
S S 
where Pmax iS the maximum value of water impact pressure 
impulse; 2t, is the duration time of impulse. 


2 
Pmax > Pmax = ky” 
T 


(8b) 
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Fig. 3. Peak pressure of monolithic plate and LWPSPS 
at different water impact velocity 
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Fig. 4. Pressure time history of water impact 
for monolithic plate and LWPSPS (C1-C6) 


If the peak impulsive pressure Pmax is observed in Eqs. (7), 
the objective now is to determine the duration time of impulsive 
2t,. The pulse pressure loading time of circular plate has been 
suggested by Chuang [4] as following estimation: 
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4r 
2t.= = 
S 
. . A C, . 
where r is the radius of circular plate; C, is the sound speed 
of air (Cy.= 340 m/s). For rectangular plates, the peak pulse 
pressure loading time can be modified as [4]: 


(9a) 


2a’ +b? 
z= C, 

It is obtained that the 2t, of monolithic panel in the present 
research is nearly 4.0 ms from Eq. (9b). This estimation is 
consisted with the numerical results in Fig. 4. Furthermore, the 
duration time of impulse 2t, is less sensitive to the variation 
of water impact velocity. As shown in Fig. 4, the water impact 
loading time of LWPSPS is longer than that of solid plate. 


Based on the numerical results, a modified estimation of 21, 
for LWPSPS is given by: 


(2.1 ~2.6)Va +b? 
C, 


It is observed that the impulsive duration time 2r, for all 
cases of LWPSPS is 4.50 ms-5.50 ms from Eq. (10). 

Moreover, the results obtained from fluid materials volume 
distribution iso-surface in initial water entry period show 
that an interesting phenomenon - “local air cushion” exists 
during the initial water impact duration (see the dashed line in 
Fig. 5). But this interesting phenomenon does not exist in initial 
water impact period of monolithic plate. This local cushion 
is mainly caused by lattice cell local deformation which has 
significant effects on the water impact pressure and dynamic 
responses of LWPSPS. Because of the existence of this /ocal 
air cushion, which plays a role as a buffer device, the peak 
value of average FSI pressure of contacted-water wet face sheet 
of LWPSPS is much lower than that of monolithic plate (see 
Fig. 3). Additionally, as shown in Fig. 4b, the high frequency 
components of hydroelastic pressure curve for LWPSPS are 
excited by the local cell deformation. 


THEORETICAL ANALYSIS: 
HYDROELASTICITY ASSESSMENT AND 
DEFORMATION CALCULATION 


(9b) 


Ts 


(10) 


2 = 


Basic vibration characteristics of monolithic 
plate and LWPSPS 


In order to evaluate the hydroelastic effects in section 
4.2, the fundamental frequency f (f,) of bending vibration 
mode for solid plate (LWPSPS) with clamped supported 


Image 
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Magnification 


boundaries is calculated firstly. The fundamental frequency 
f of bending vibration modes of a clamped monolithic plate 
[30] is given by: 


where m, is the mass unit area, D is the bending stiffness of 
plate, a and b are the length and the width of plate. In the present 
analysis, the fundamental frequency f of this monolithic plate 
is 938 Hz and the analytical result is in good agreement with 
the FEM simulation result (922 Hz). 

Because of the discrete characteristics of lattice truss 
core, the calculation of the natural frequency of LWPSPS is 
much more complicated. Usually, the generalized equivalent 
model can be mostly used in the preliminary analysis stage of 
design process. It can be used to reduce the time spent for the 
analysis of LWPSPS. The homogenization method also makes 
it possible to avoid the problems involved in heterogeneities. 
Here, a new equivalent method to calculate the fundamental 
bending vibration frequencies of LWPSPS is proposed as 
follows (Fig. 6): 

- Step 1: calculate the equivalent shear rigidity of LWPSPS 
lattice core. 

There are a few methods to calculate the equivalent shear 

constant of lattice materials. Liu [31] proposed a matrix 

displacement method of lattice truss based on the small 
displacement assumption. Therefore, this homogenization 
method is employed here to calculate the shear rigidity of 

LWPSPS lattice core as follows: 


2t? 

E Ea (12a) 
ree I’ cos’@ sing 

E, = EÐ, sinf ọ (12b) 

2 
r | H. +(t,+t,)/2 ] o E io 
2E  — J c 

° H, Cis C4 


where p, is the relative density of this lattice core [32], 
E, is the normal modulus of LWPSPS core [32], Cois 
the transverse shear rigidity matrix of LWPSPS core (the 
detailed expressions of the homogenized moduli of lattice 
core Ch, , (a, B =4, 5) can be found in literature [31]), the 
definitions of 1, ọ and t, can be seen in Fig. 1 and p, is the 
material density of the cell structure. 
- Step 2: calculate the bending rigidity of LWPSPS. 

For the thin face sheets of LWPSPS, only bending 
deformation is considered (the bending stiffness D, of face 


Image Magnification 


Local air cushion 
Fig. 5. Fluid volume distribution of initial entry period of LWPSPS and free surface contour iso-surface 


Fig. 6. Equivalent model of LWPSPS 
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Fig. 7. The 1st - 3st bending vibration modals of LWPSPS (Case1) 


sheets is included). The bending stiffness of sandwich panel 
(LWPSPS) with face layers of non-equal thickness [33] is 
given by: 


E((H, +t, +t,) -H’) 


1 2 
2 
i a(t; =f.) 
(H. +t, +t.) 
H, ct lg Thy (t, +t,) 


+— + 
(H, +t, +t,) H, 


The meanings of all nomenclatures are given in Fig. 1 and 

Fig. 6. 

- Step 3: calculate the natural frequency of LWPSPS with 
clamped supported boundary conditions. 
If the bending rigidity and shear rigidity of LWPSPS are 
obtained in step | and step 2, the last step is to calculate 
the natural frequency with clamped supported boundary 
conditions. In this key step, an approximate method [34] 
is adopted to analyze the fundamental vibration frequency 
of the LWPSPS. Due to the complexity of approximate 
mathematic equations, only the calculation results are given 
here. Detailed steps of solving approximate equations can 
be found in literature [34-35]. 


The first modal frequency f, of LWPSPS (Case 1) based on 
the equivalent-approximate method is 1698 Hz which is close 
to the FEM result (1553 Hz) in Fig. 7 (the error is 9.34 %). It 
must be pointed out that the equivalent-approximate method is 
more accurate in high frequencies [35]. Thus, the error between 
FEM numerical result and equivalent approximate result for the 
second (third) bending vibration frequency is less than 5 %. 


Analysis of hydroelastic effect 
As mentioned in the previous section, the transient response 


of a free fall structure subjected to hull-water impact loading is 
greatly complicated. The evaluation of hydroelastic (FSD effect 


is an important part of this investigation topic. In Ref. [8], the 
hydroelastic effect indicator is formulated by Bereznitski: 


Rg = Typp/Typp (14) 


where T, pg is the water impact loading time, and T,,, is the first 
dry natural period of vibration. Compared Eq. (14) with the 
definition 2q, f, (or 2t, f), it is found that these two expressions 
have the same physical meaning. Bereznitski [8] concluded 
that the hydroelastic effects become important when a ratio Rg 
(or 21, f.) < Rpm Where Rpiim ~ 2.0. 

Jones also pointed out that the responses of structure can 
be seen as quasi-static if the duration of water impact loading 
time 2t, is much longer than the first natural bending vibration 
period T (2t, >> T) [27]. In fact, for a single-degree spring- 
mass system (SDOF), the deformations calculated for the same 
load applied statically agree with the corresponding maximum 
dynamic deformations to within approximately 17 % when 
2t/T > 1.75 [36]. 

According to Eq. (11), the responses of this clamped plate 
can be seen as quasi-static (21, f= 4.2 - 5.2 > 2.0). For LWPSPS, 
the multiplication 21,f, (21, = 4.50 - 5.50 ms, f, = 1560 Hz) is 
6.8-8.3. Thus, the dynamic responses of LWPSPS in the present 
case can also be seen as quasi-static. Normally, this criterion is 
absolutely true as a matter of fact that the slamming duration 
time is about 0.02 s - 0.1 s in practical case [19], which is much 
longer than that of the present investigation. 

While the dynamic inertial effect cannot be ignored, the 
calculation of dynamic deflection in linear elastic range under 
uniform impact pressure is given by following equation [36]: 


K(t)= wa aa [p@sinott —t)dt (15) 


where K(t) is the dynamic inertial effect factor function, @ is 
the first circular frequency of structure, Pmax 1S the peak pressure 
value of impact, p(t) is the formula of impact loading, w,(t) is 
the load-deflection function and w, is the static deflection of 
structure under uniform maximum pressure Pax The maximum 
value of K(t) at t = tẹ is commonly considered in the practical 
case, which satisfies the following conditions: 
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dK (t) _ d°K(t) <0 
2 t= 


dt hev» dt ~ 


In this study, the maximum value of K(t) (Here it is defined 
as 1,4) can be easily calculated for monolithic plate based 
on Eqs. (7)-(10) and Eqs. (15)-(16). The value of K is only 
associated with the fundamental circular frequency œ and 
the impact duration time 2t, in Eq. (15). If the fundamental 
frequency is obtained in the previous section 4.1, the maximum 
value of K(t) for LWPSPS (Here it is defined as n) can be 
calculated from Eqs. (15)-(16). 

In the present investigation, n,a (for monolithic plate) is 
equal to 1.20 - 1.25. For LWPSPS, n,a is equal to 1.05 - 1.10. 
Compared with computation results, the values for n,a and 
Na are both verified reasonable. Therefore, the values of n,a 
(1.20 ~ 1.25) and n, (1.05 ~ 1.10) indicate that the dynamic 
inertial effects is not remarkable during the water impact 
process in these cases. 

Considering the added mass of the first vibration modal of 
structure for a v-shaped ship section, a similar definition given 
by Stenius et al. [6, 10, 37] is defined as follow: 


R=T,p/Typ (17) 


where T, pis the loading period based on the Wagner theory, and 
Typ is the first wet natural period of vibration. It is concluded 
that the hydroelastic effects are negligible and the response can 
be set as quasi-static for R > 4 [6, 10, 37]. 


Obviously, the first wet natural period of vibration Typ 
is longer than the first dry natural period of vibration Typp,. 
Though a detailed hydroelastic analysis of two-dimensional 
v-shaped transverse section was conducted in [6, 10, 37], the 
hydroelastic problem of a flat panel with zero wedge deadrise 
angle are neglected in previous investigations. 

The indicator of hydroelastic effect introduced by 
Bereznitski [8] and the inertia effect indicator introduced by 
Jones [27] are both based on the dry natural period of vibration. 
However, the added water mass is not considered in their 
analysis. Pointed out by Stenius et al. [6, 10, 37], hydroelastic 
effects incorporate both dynamic inertial effect and added mass 
effect. Thus, the added water mass should be considered in the 
evaluation of hydroelastic effects. 

For a rectangular clamped monolithic plate with one face 
contacted with water, the added water mass coefficient is 
defined as f,/f according to Ref. [30]: 

fa 1 


1+6 = ae) 
yt 

where f is the first dry natural frequency of vibration, f, is the 
first wet natural frequency of vibration, t is the thickness of this 
rectangular monolithic panel, y, is the specific gravity of steel 
(here y,.= 7.85), a is the length of this monolithic plate and ¢ 
is the constant based on the length-width ratio of a/b. In the 
present case, the ratio of a/b is equal to 2.2 and the value f, /fis 
equal to 1.34. According to this observation, the first wet natural 
frequency of vibration f, is equal to 699 Hz which is consistent 
with the result (677 Hz) obtained by FSI vibration analysis in 
ANSYS by using acoustic fluid element Fluid30 [38]. 


As the FSI occurs, the calculation of wet natural frequencies 
of LWPSPS (here it is defined as f,s) is greatly complicated. 
Up to now, no study analyzed the FSI vibration characteristic 
of LWPSPS. Thus, the acoustic-solid coupling method in 
commercial code ANSYS is adopted to calculate the wet 
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natural frequencies of LWPSPS. The fluid is also modeled 
as non-viscous acoustic fluid element Fluid30, and Lagrange 
elements (shell 63 and beam 181) are adopted for solid parts 
(face sheets and core trusses); where the sum of all elements is 
nearly 80,000. When the added water mass effect is taken into 
account, the first wet natural frequency of LWPSPS is 1034 Hz 
for Case 1 - Case 6, which is about 2/3 of the first dry natural 
frequency (1553 Hz). 

Once the first wet natural frequency of LWPSPS is obtained, 
the hydroelastic indicator of R can be calculated. By applying 
Eq. (14), the value of hydroelastic indicator R for LWPSPS 
(C1-C6) is 4.653 - 5.687 (2t, = 4.50 - 5.50 ms). In Ref. [6, 
10, 37], Stenius et al. concluded that the hydroelastic effects 
are negligible when R > 4. Thus, one could give a definitive 
conclusion that the hydroelastic effect of LWPSPS is not 
significant based on this criteria. 

This conclusion can be validated from 3D FEM FSI 
simulation results (3D-FEM-FSI), 3D FEM dynamics 
simulation results (the FSI effects are not considered, 3D- 
FEM-Dynamic) and 3D FEM quasi-static simulation results 
(3D-FEM-Static) in Fig. 8. The comparison of these analyses 
shows that the results of 3D-FEM-FSI and 3D-FEM-Dynamic 
are very close, which means the hydroelastic effect is not 
significant. As shown in Fig. 8, due to the neglect of inertial 
effect, the quasi-static results are slightly lower than the results 
of dynamic calculation. 
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Fig. 8. Comparison of 3D-FEM-FSI, 3D-FEM-Dynamic 
and 3D-FEM-Static Results 


4.3. Engineering approach of maximum deformation for 
LWPSPS 


Next, the analytical investigations are performed to study 
the deformation characteristics of LWPSPS. As shown in 
Fig. 9, the inflection points exist at adjacent lattice core cells. 
The total deformation field of LWPSPS can be divided into 
two components: 

1) The global bending deformation field of LWPSPS 


W max.T)- 
2) The local bending deformation field of each cell (Wmax 1)- 


Based on this deformation field assumption, a mathematical 
model to describe the total maximum deformation field is 
presented as follow: 


Ww Ww + Ww 


(19) 


max max,T max,L 


- Calculation of the global bending deformation W maxr 
The global linear elastic maximum deformation w,,,,7 of 
sandwich plate with clamped boundary condition under 
static uniform pressure is given by [35]: 
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where C, is the equivalent shear stiffness of lattice core and 


Pmax 1S the maximum value of average FSI pressure, the other 
nomenclatures have the same meanings as mentioned in the 
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Fig. 9. The out-plane deflection counters plot of LWPSPS 


previous section. Firstly, the bending and shear stiffness of 
lattice core D, and C, are calculated from Eqs. (12) - (13). 
Then the coefficient H,, and G, are obtained from Eq. (20b) 
and Eq. (20c). If all the parameters in Eqs. (20b) - (20e) are 
observed, the maximum static deformation can be easily 
solved from Eq. (20a). 
Calculation of the local bending deformation field of each 
cell Wnai 
The confirmation of boundary condition is a key step to 
analyze the local deformation field of each cell. Here, three 
boundary condition assumptions of lattice cell are compared 
as follows: 
1. Cylindrical bending with built-in edges (Analytical-C, see 
Fig. 10a)[39] 


4 
Pmnax“ c 
W. -e 21 
Da 384D, (eaa 


N 


Cylindrical bending with simply supported edges 
(Analytical-S, see Fig. 10b) [39] 


— 5 Prado 


W. = 
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3. Bending ofeach lattice cell supported by rows ofequidistant 
trusses (Analytical-T, see Fig. 10c) [39] 


(21b) 
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Where w,,,x; 18 the local maximum deformation, Pmax is the 
maximum value of water impact pressure loading on central 
cells, d, is the width of cell (see Fig. 1) and D, is the flexural 
rigidity of cell face sheet defined as: 


O Eg 
12 -v 


If the inertial effect during water impact loading process 
is significant (2t, f, ~ 1.0), the total static deformation W max.r 
of LWPSPS must be multiplied the dynamic load inertial 
coefficient 1,,°. The value of n, is equal to 1.05 - 1.10 for 
case | - case 6 from Eq. (16) and Eq. (17). Here, the value 1.10 
of n, is adopted in the following investigation. 

Compared with the 3D-FEM FSI simulation, the maximum 
deformations Wax of LWPSPS under different water impact 
velocity for the three assumptions (Analytical-C, Analytical-S 
and Analytical-T) are illustrated in Fig. 11. It is clear that the 
analytical deflections match well with computed results by 


D, (22) 
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Bottom jes sheet 


b) Analytical-S 


Y 


Bottom face sheet 


c) Analytical-T 
Fig. 10. Three local cell deformation field assumptions of LWPSPS 


adopting the Analytical-S assumption. The results show that 
the maximum percentage difference (100% x(Wsim - Want)/Wsim) 
between this engineering approach (W=!) and the FSI simulation 
results (Ws™) is less than 9 % (impact velocity v = 3 - 6 m/s). 
Furthermore, the results based on Analytical-T assumption 
can give a better estimation in the low velocity range (impact 
velocity v = | - 2 m/s). The maximum difference between the 
simulation and the approximate method is less than 3 %. As 
shown in Fig. 11, the results based on Analytical-C assumption 
for LWPSPS are unacceptable. 
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Fig. 11. Comparison of the maximum deformation of LWPSPS by using 
approximate method and FSI simulation 


PARAMETRIC STUDY 


In this section, the effects of design parameters of t,, ts 
H. txt, on water impact pressure and maximum deflection 
of LWPSPS are examined by numerical simulation and 
engineering prediction. 


Effects of thickness 
of face sheets t, and t, (C7-C14) 


To investigate the effects of thickness of face sheet on water 
impact loading pressure, LWPSPS with four different bottom 
face sheet thicknesses (C7 - C10) and four different top face 
sheet thicknesses (C11 - C14) are considered. The material 
properties, plate geometry, water impact velocity and boundary 
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conditions are the same as those that adopted in the previous 
study for Cl - C6. 

With respect to the fundamental case C1, the variation in 
thickness of bottom face sheet t, leads to the significant shifts 
of water impact loading pressure and maximum deflection, 
as shown in Fig. 12 and Fig. 13. However, the influence of 
variation in thickness of top faces t, is insignificant. With 
the increase in thickness of t, (the contact water face sheet 
“hardens’”’), the FSI peak pressure of central point increases, 
as illustrated in Fig. 12. On the contrary, the average FSI peak 
pressure of contacted water face sheet has a minimum when t, 
is nearly equal to 4.5 mm. 
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Fig. 12. Peak water impact pressure magnitudes at central point and 
average FSI area with different face sheet thickness: top and bottom 
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Fig. 13. Maximum deformations of LWPSPS with different bottom face sheet 
thickness predicted by analytical method and FEM 


It can be seen clearly from Fig. 13 that, with the thickness of 
bottom face sheet increases, the deformation curves demonstrate 
a progressive decrease when t, < 4.5 mm. Contrarily, when 
t, > 4.5 mm, a slowly increasing trend of deformation curve 
is observed. In an interpreting way, when t, > 4.5 mm, though 
the flexural rigidity of LWPSPS increase slowly, both the 
central pressure and the average pressure increase faster, as 
shown in Fig. 12. 

With no remarkable change of water impact pressure, the 
maximum deformation magnitude of LWPSPS is smaller when 
t, becomes thicker (Fig. 14). The approximate method also gives 
good estimations here under the Analytical-S assumption, and 
the error is less than 10% in most of the cases. With very few 
exceptions (t,or t, = 2.219 mm), the precision of this analytical 
prediction is high enough. As a possible factor, the Analytical-S 
assumption cannot reflect the true boundary condition of local 
lattice cell when the face sheet is relatively thin. 
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Fig. 14. Maximum deformations of LWPSPS with different top face sheet 
thickness predicted by analytical method and FEM 


Effects of core thickness 
of LWPSPS H,(C15-C18) 


Next, the influence of core thickness on the water impact 
pressure and maximum deflection of LWPSPS is considered 
(C15 - C18, see Table 1). The peak pressure magnitude versus 
core thickness for C15 - C18 under the same water impact 
velocity 5m/s are demonstrated in Fig. 15. Along with the 
increasing thickness of core, the peak pressure value moves 
to a lower magnitude. As illustrated in Fig. 15, the slope 
of decreasing curve changes slowly for central pressure 
when H, is in the range of 1.825 - 2.825 cm. But the slope 
of central pressure curve drops dramatically in the range of 
H, = 2.825 - 5.825cm. Moreover, a slight decrease of 
average FSI peak pressure magnitude with a thicker core is 
observed. 

To quantify the effects of core thickness of LWPSPS on 
the structural response, the predicted deflections versus core 
thickness are plotted for these four cases with three different 
local boundary condition assumptions, as shown in Fig. 16. As 
the thickness of core increases, it is shown that the deformation 
of LWPSPS decreases rapidly in the range of H, = 1.825 - 
4.825 cm. This is because of the decrease of FSI pressure and 
the increase of bending rigidity of LWPSPS with thicker core. 
The best prediction of maximum deformation of LWPSPS 
is also obtained under the Analytical-S assumption. When 
H, > 2.5cm, the engineering approach results are consistent well 
with those obtained from numerical calculation. Also with very 
few exceptions (when H, < 2.5cm), the discrepancy between 
the numerical results and the analytical solutions displays in 


Fig. 16. Compared with FEM result, the error of the prediction 
is about 25 %. 
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Fig. 15. Peak water impact pressure magnitudes at central point and 
average FSI area with different core thicknesses 
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Fig. 16. Maximum deformations of LWPSPS with different core thicknesses 
predicted by analytical method and FEM 


Effects of cross section area t,xt, (C19 - C22) 


Finally, the influence of the section area of lattice column 
as an indicator of relative density is discussed. According to the 
calculation results, nearly 2 times increasing of relative density 
of lattice core (from 4.83 to 9.47 %) results in 10 % decrement 
of FSI pressure located in the central point and 10 % increment 
of the average FSI pressure. 

As demonstrated previously, the contacted water face 
sheet and core thickness both have significant effects on the 
dynamic response behavior of LWPSPS. But it is shown 
that the change of relative density has a slight effect on the 
dynamic response behavior. In contrast, in the evaluation of 
the energy absorption characteristic of LWPSPS under intense 
impact loading such as UNDEX (underwater explosion) or 
AIREX (air explosion), the relatively density of lattice core 
plays a pivotal role on the dynamic mechanical properties 
of lattice materials [17-18]. However, unlike with the plastic 
collapse under these intense impact loads (UNDEX or 
AIREX), the serious failure modes (like the plastic dynamic 
buckling) would not occur under the practical water impact 
loading (here, impact velocity v = 5 m/s) only if the relative 
density is extremely low. 

The maximum deformations of C19 - C22 are also 
predicted by numerical simulation and analytical method. 
With no remarkable change of the water impact pressure, 
the corresponding results demonstrate that the maximum 
deformation decreases slowly (about 20 %). 
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CONCLUSIONS 


The hydroelastic characteristics for light weight pyramidal 
sandwich plate structures (LWPSPS) under water impact by 
using 3D simulation are studied. A theoretical hydroelasticity 
assessment model for LWPSPS and a new engineering approach 
method to predict the maximum deformation of LWPSPS under 
water impact load are developed. 

In numerical analysis, an arbitrary Lagrangian and Eulerian 
algorithm (ALE) in the nonlinear finite element code LS-DYNA 
is used to simulate the fluid-solid interaction phenomenon 
between structure hull and fluid medium. The equivalent 
core model and the approximation method are adopted to 
calculate the fundamental bending vibration frequencies of 
LWPSPS. Furthermore, based on the first fundamental modal 
frequency, the hydroelastic characteristic is discussed in detail 
by introducing a non-dimensional criterion 21,f,. By dividing 
the real displacement field of LWPSPS into global displacement 
field and local displacement field, the real maximum 
deformation solutions are obtained. Compared with 3D FEM 
results, the results of the engineering calculation method 
presented in this analysis are obtained with good accuracy. 
At last, the influences of key design parameters of LWPSPS 
on dynamic behavior are quantified, including the face sheet 
thickness, the thickness of lattice core and the relative density 
of LWPSPS core. 

The following conclusions from the viewpoint of 
hydroelastic characteristics of LWPSPS can be drawn: 

1) In the form of water impact velocity, the performance 
of water impact pressure for LWPSPS is close to two 
parameters parabolic mathematical expression. Compared 
with the monolithic plate, in medium-high water impact 
velocity range (2m/s-6m/s), the advantage of deceasing 
water impact pressure peak for LWPSPS is significant. A 
notable phenomenon called “local air cushion” exists during 
the initial water impact process of LWPSPS caused by the 
cell local deformation, which has significant effects on the 
dynamic characteristic of LWPSPS. 

2) The hydroelastic effect of LWPSPS is not significant in the 
present study because the non-dimensional criteria 2t f>4, 
which means the duration time of water impact pressure 
loading is much longer than the first wet natural period of 
vibration. 

3) Results from the proposed engineering analytical model 
are in good agreement with those from finite element 
predictions for LWPSPS. All the geometric parameters are 
included and the local structural displacement is captured 
in this model. The predictions based on this method show 
that, compared with the other two local deformation field 
assumptions (cylindrical bending with simply supported 
edges and cylindrical bending with built-in edges), the 
local cylindrical bending assumption with simply supported 
edges each row lattice core is more realistic to reflect the 
true boundary condition in most cases. 

4) In the parametric study of LWPSPS, it is found that the 
thinning of contacted water face sheet reduces the central 
water impact loading pressure, while the average FSI water 
impact loading pressure has a local minimum value at 
a certain thickness. Additionally, along with the increasing 
of core thickness, the peak pressure moves to a lower 
magnitude. Another key design parameter-the relative 
density of lattice core, has no significant effect on the FSI 
pressure and maximum deflection when the relative density 
of core is doubled. 
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Graphical interpretation of the power of energy 
losses and power developed in the hydrostatic 
drive and control system elements 
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ABSTRACT 


This publication presents, analyses and compares the areas of the power fields of energy 
losses occurring in the elements of two hydraulic systems with different structures of the 
hydraulic linear motor speed control. 


Keywords: hydraulic motor; energy balance; energy losses; energy efficiency; 
capacity pump; hydraulic servo mechanism; hydrostatic drive system 


INTRODUCTION 


The paper presents graphical interpretation of the power 
of energy losses and power developed in the elements of 
two different systems with the throttling series control of the 
hydraulic linear motor speed. The analysis was carried out 
by comparing, at selected hydraulic linear motor operating 
parameters, the areas of the power fields of energy losses 
occurring in the elements of those structures. 

The investigations involved two systems with proportional 
directional valve, fed by a constant capacity pump (Fig. 1): 
a) with the use of an overflow valve - a constant pressure 

structure, 

b) with the use ofa hydraulic cylinder supply conduit pressure 
controlled overflow valve - a variable pressure structure. 


The investigated structures operated at the same hydraulic 
linear motor operating parameters, i.e. its load F,, and 
speed Vy- 

The analysis allows to compare the values of power AP 
of losses ensuing from the used structure of control of the 
hydraulic linear motor speed as well as the value of power Pp. 
absorbed (consumed) by the pump from its driving electric 
motor, power necessary for providing the required stable value 
of useful power Pu, = Fy Vy of the hydraulic linear motor 
driven by the pump. 

Graphical presentation, by fields of specified areas, of the 
power of energy losses in the hydrostatic drive and control 
system elements as well as power developed in the hydraulic 
displacement machines used in the system, becomes a tool 
facilitating comparing the values of particular losses [1]. 
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Fig. 1. System with proportional directional valve fed by a constant capacity 
pump [3]: a) with the use of an overflow valve — p = const structure, 
b) with the use of a hydraulic cylinder supply conduit pressure 
controlled overflow valve — p = var structure 


Presentation of the fields of power of energy losses allows to 
make conclusions regarding e.g. elimination of the power of 
structural volumetric and pressure losses in the motor speed 
throttling control elements, in the proportional control systems 
and in the hydraulic servo-mechanism systems. Graphical 
interpretation by the field areas of the power of energy losses 
in the hydrostatic drive system elements and of the power 
developed by the system elements allows to compare those 
losses and powers with the area of the reference power field 
defined by the product Qp, p, of the theoretical pump capacity 
and the system nominal pressure [1]. 


HYDRAULIC SYSTEM WITH 
PROPORTIONAL CONTROL OF THE 
CYLINDER FED BY A CONSTANT 
CAPACITY PUMP IN THE CONSTANT 
(P = CONST) AND VARIABLE (P = VAR) 
PRESSURE SYSTEM 


Proportional control of a cylinder consists in throttling the 
liquid stream both at its inlet and outlet [3]. 

The basic proportional control system is a system fed by 
the constant capacity pump. The overflow valve SP (Fig. la) 
determines the system nominal pressure. The pressure decrease 
in the cylinder compensates the load on the cylinder. The 
proportional directional valve generates two pressure drops 
at the cylinder inlet and outlet. The pump in the p = const 
system must generate, before the overflow valve, pressure 
not lower than pressure required by the cylinder. Therefore, 
the hydraulic cylinder or the system working cylinder may 
require pressure, depending on the load, in the range from zero 
to the nominal value. When the load approaches the nominal 
value, pressure decrease in the directional valve throttling slots 
tends to zero. It may be said that the pump-overflow valve 
assembly in the p = const system is ready to feed the system 
with the maximum pressure and maximum capacity, but most 
often it is not used to that extent as the working element is 
loaded with a force that requires pressure drop smaller than the 


Q 


OQ = pi Meo 


o| 


nominal value [3]. A constant pressure system achieves a high 
energy efficiency,equal to the efficiency of a system without 
throttling control, only at the point of maximum values of the 
controlled hydraulic linear motor load coefficient M,, and 
speed coefficient Oy. The system efficiency n decreases rapidly 
with decreasing motor load and particularly with simultaneous 
decreasing motor speed [2]. 

The variable pressure (p = var) structure is represented 
by a system with constant capacity pump cooperating with 
an overflow valve controlled by the cylinder inlet pressure 
(Fig. 1b). This is an advantageous solution from the viewpoint 
of the cylinder energy efficiency as well as of the pump and 
the whole control system efficiency. The variable pressure 
(p = var) structure with the overflow valve SPS controlled 
by the current directional valve outflow to cylinder pressure 
allows to adjust the pump discharge conduit pressure to the 
current cylinder load, which limits the pressure loss in the 
working liquid outflow slot from the directional valve to 
the tank. Additionally, the system maintains constant piston 
speed irrespective of the load. This is an effect of maintaining 
practically constant pressure drop App,, in the proportional 
directional valve throttling slot [3]. 

Fig. 2 presents graphical interpretation of the power 
of energy losses in elements of an individual system with 
proportional control of a hydraulic cylinder, fed by a constant 
capacity pump cooperating with an overflow valve in a constant 
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Fig. 2. Graphical interpretation of the power of losses in the hydrostatic drive and control system elements. An individual system with series throttling 
control of the hydraulic linear motor speed, fed by a constant capacity pump cooperating with the overflow valve in a constant pressure system — p = const; 
the series throttling control assembly in the form of [3]: 

- adjustable throttling valve; - adjustable two-way flow regulator; - servovalve; - proportional directional valve 
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pressure system (p = const), and Fig. 3 — with an overflow valve 

controlled in a variable pressure system (p = var). 

The pump operation nominal pressure p, is determined 
by the need to ensure the hydraulic linear motor a maximum 
pressure drop Apyimax to cope with the maximum force F ymax 
on the cylinder piston rod that may occur from time to time as 
an effect of the cylinder load. 

The current cylinder useful power Puu = Fy Vy is a product 
of the cylinder piston rod force F,, and the piston rod speed. 
The hydraulic linear motor useful power P,,, depends on the 
current load and is independent of the control structure and of 
the losses in the elements of a hydrostatic drive system with 
a specific structure [1]. 

In figure 2, the hydraulic linear motor current useful power 
Pu = Fy Vy 1s presented as the area of the white rectangle, to 
which the following fields are “added”: 

- field APum = Fum Vm Of the power of mechanical losses in 
the cylinder; 

- field AP w = APyi Qu, Of the power of volumetric losses in 
the cylinder, 

- field AP\,, = APmp Qm of the power of pressure losses in the 
cylinder, 

- field AP. = Apc Qu of the power of pressure losses in the 
system conduits, 

- field AP,,, = App: Qu of the power of structural pressure 
losses in the throttling control assembly (in the proportional 
directional valve), 

- field AP... = psp (Qp - Qu) of the power of structural 
volumetric losses in the throttling control assembly (in the 
overflow valve), 

- field AP,, = Apy, Qu of the power of pressure losses in the 


- field AP,, = APppi Qp, of the power of volumetric losses in 
the pump, 

- field AP,,, = Mpm ©p of the power of mechanical losses in 
the pump. 


The sum of areas of the cylinder current useful power 
Pu rectangle and the AP rectangles representing the values 
of power of particular losses occurring at a given instant in 
the hydrostatic drive and control system elements, makes up 
the area of rectangle corresponding to the current power Pp. 
absorbed (consumed) by the pump from its driving electric 
motor and equal to the product of the current torque M, and 
current pump shaft angular speed — Pp. = M,@, [1]. 

The power P,, absorbed by the pump from the driving motor 
may be greater than the reference power p, Q,,—a product of the 
nominal pressure p, and pump theoretical capacity Qp [1]. 

In figures 2 and 3 the pump capacity is represented by 
two descending lines originating at point Q,,. The higher line 
illustrates the pump capacity in the situation when the pump 
volumetric loss coefficient k, has a zero value (k, = 0). The line 
below represents the pump capacity at k, > 0. The volumetric 
losses occur in the latter situation [3]. 

The Q,, curve presents the proportional directional valve 
characteristic with a given pressure drop App, in the directional 
valve and with a given throttling slot area fpg (Fig. 2 and 3). 
At point “A” it crosses the overflow valve (SP) characteristic 
Psr = f(Qp — Qu). In effect, the intensity Qu of flow through 
the throttling slot to the cylinder is obtained, and with a given 
piston and piston rod area - the speed v,,. The working point 
A is a result of the overflow valve SP and of the proportional 
directional valve characteristics [3]. 

The cylinder useful power Pw. is a product of its speed 
Vu and load Fy. In other words, it is power developed by the 
cylinder on the piston rod. The cylinder useful power Py, field 
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is marked white to distinguish it from the power of losses in 
the system. 

Power AP u, of volumetric losses is a function of the pressure 
drop Apy; in the cylinder. 

Power AP,,,, of mechanical losses is a function of the load 
force Fu- 

Power AP, of losses in the system conduits is a product 
of the sum of resistance of flow Ap, and flow intensity Qu 
towards the cylinder. 

Power AP.,, of structural pressure losses is a product of the 
sum of pressure losses App, in the directional valve throttling 
slots and of the flow intensity Q,, corresponding to the cylinder 
speed v,,. It can be reduced almost to zero during cylinder 
operation at the F,,,,.. load. 

Power AP,,, of structural volumetric losses is a product of 
pressure psp in the pump discharge conduit and of the flow 
intensity Q, directed, through the overflow valve SP or the 
controlled overflow valve SPS, to the tank [3]. It decreases 
almost to zero when the cylinder operates with maximum 
speed Vyuma 

There are three types of losses in the pump. 

Power AP,,, of mechanical losses in the pump occurs 
between the pump working chamber and the shaft. It is 
proportional to the pump capacity per one revolution qp and 
to the pressure increase Ap; in its working chambers. The App; 
value is influenced by the resistance of flow App, in the pump 
channels [3]. Fig. 2 presents the power of pump mechanical 
losses in the p = const system as a field whose width is 
determined by the torque M,,, of mechanical losses and height 
corresponds to the pump shaft angular speed œp. 

In the p = var system (Fig. 1b and 3), the pump working 
pressure p,,, controlled by the controlled overflow valve SPS, 
is set at a level higher by the value: 


APsps = Pl me APcimax=Cte 


than the current pressure p, in the throttling directional valve 
outlet conduit to the hydraulic linear motor. The value Apsps of 
the pressure difference APsps = Pp. — po Must ensure obtaining, 
with the throttling directional valve slot DE1 controlling the 
flow intensity Q,, feeding the hydraulic linear motor, the flow 
intensity Qu equal to the theoretical pump capacity Qp, Le. 
Qu = Qa. The slot cross-section area DE1 reaches then its 
maximum value fpgimax » With a possibility of obtaining the 
AP pg I|fpe1max Qp, Op required by the throttling directional valve 
structure and at the same time the capability of overcoming the 
maximum resistance of flow APcima that may appear between 
the pump and the directional valve. The pressure value p, 
before the throttling directional valve slot DE1 is equal to p, 
= Pr — Ape, [1]. 

The current value of the pump discharge pressure pp, higher 
by Apgps than the current p, value at the throttling directional 
valve outlet to the hydraulic linear motor, results from the 
pressure value py required by the hydraulic cylinder at its inlet. 
The maximum limit pressure value Ppzmax 1n the pump discharge 
conduit is determined by the overflow valve SP, whose opening 
pressure Pspo is equal to the system nominal pressure p, [1]. 

In the p = var system (Fig. 3), the width of the field of power 
AP,,, of mechanical losses in the pump is proportional to the 
pressure and is smaller as the pressure in the pump working 
chambers is lower. In effect, it requires a smaller torque on the 
pump driving motor shaft [3]. 

The field AP,, of the power of pressure losses in the p = var 
system pump is greater than the p = const system AP,,. 

The power AP,, of volumetric losses in a p = var system 
pump is distinctly smaller than the power AP,, of those losses 
in a p = const system pump. 
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Fig. 3. Graphical interpretation of the power of losses in the hydrostatic drive and control system elements. An individual system with series throttling 
control of the hydraulic linear motor speed, fed by a constant capacity pump cooperating with the overflow valve in a variable pressure system — p = var; the 
series throttling control assembly in the form of [3]: 

- adjustable throttling valve (together with a controlled overflow valve, creating a three-way flow regulator); - servovalve; - proportional directional valve 


The intensity Qu of flow towards the hydraulic cylinder 
results from the pressure drop Ap,,, and Ap,,, in the 
proportional directional valve slots. On the other hand, those 
pressure drops change when the pressure drop Ap,, required by 
the hydraulic cylinder changes [3]. 


COMPARISON OF THE AREAS OF FIELDS 
OF THE POWER OF ENERGY LOSSES IN 
THE ELEMENTS OF TWO STRUCTURES 

— P= CONST AND P= VAR 


Figures 4 to 9 present graphical interpretation of the power 
of losses in the hydrostatic drive and control system elements, 
investigated in a constant pressure (p = const) and variable 
pressure (p = var) system, with the hydraulic cylinder load 
values Fu = 10 KN, 25 KN and speed values v,, = 0.05 m/s, 
0.20 m/s and 0.35 m/s. The fields of power of losses AP in the 
constant pressure system are shaded red and fields of power of 
losses in the variable pressure system - shaded blue. 

At the cylinder zero load F,, = 0 there is no field of the 
cylinder useful power Py, and only losses occur in the p = const 
and p = var system. 

In Fig. 4, at the cylinder load F,,= 10 kN and speed 
Vy = 0.05 m/s, great structural volumetric losses AP,,, are 
visible in the SP overflow valve, as with small cylinder speed 


the flow intensity Qu towards the cylinder is small. Most part 
of the flow intensity Q, generated by the pump is directed to 
the tank as the Q, intensity. 

Fig. 5 illustrating a variable pressure (p = var) structure 
shows structural volumetric losses AP,,,, occurring in the 
controlled overflow valve SPS, the greatest among the losses 
as with small cylinder speed the flow intensity Qu towards the 
cylinder is small. Most part of the flow intensity Q, generated 
by the pump is directed to the tank in the form ofa Q, intensity. 
At the same time, power AP „ of volumetric losses is much 
smaller here than AP,,, in the p = const system (Fig. 4), as the 
pump working pressure pp, set with the controlled overflow 
valve, is much lower. 

Comparing the energy balance of the p = var system shown 
in Fig. 5 with the p = const system balance shown in Fig. 4, 
with the same cylinder load and speed coefficients, one can 
notice much smaller structural losses AP,, in the proportional 
directional valve and in the controlled overflow valve as well as 
the pump volumetric losses AP,, in the p = var system. This is 
connected with adjusting the p,, pressure in the pump discharge 
conduit to the current cylinder load, i.e. to the p, pressure. As 
Fig. 5 illustrates the energy balance of a p = var system under 
a small load, losses in the system are smaller. 

Losses in the conduits are the smallest in this case, as the 
power of losses in the conduits AP, is a product of the sum 
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Fig. 4. Energy balance of a constant pressure (p = const) system investigated at the cylinder load FM = 10kN and speed vM = 0.05m/s 
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Fig. 5. Energy balance of a variable pressure (p = var) system investigated at the cylinder load FM = 10KN and speed vM = 0.05m/s 
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Fig. 6. Energy balance of a constant pressure (p = const) system investigated at the cylinder load FM = 25kN and speed vM = 0.20m/s 
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Fig. 7. Energy balance of a variable pressure (p = var) system investigated at the cylinder load FM = 25kN and speed vM = 0.20m/s 
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Fig. 8. Energy balance of a constant pressure (p = const) system investigated at the cylinder load FM 
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Fig. 9. Energy balance of a variable pressure (p = var) system investigated at the cylinder load FM 
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of resistance of flow Ap, and flow intensity Qu towards the 
hydraulic cylinder, and that intensity is small due to small 
cylinder speed v,, = 0.05 m/s. 

Figures 6 and 7 present the field P,,, of the hydraulic 
cylinder useful power (white field) and the fields of the power 
AP of losses that occur in the p = const system (Fig. 6) and 
p = var system (Fig. 7) at the cylinder speed v,, = 0.20 m/s and 
cylinder load F,,= 25 kN. The fields of power AP of losses are 
shaded red. 

The share of structural volumetric losses AP „ and pressure 
losses AP,,,, in the total energy balance of a constant pressure 
system is relatively great compared with other losses in the 
elements. The pressure energy accumulated in the liquid 
pressed by the pump is only partially used effectively in 
the cylinder. The remaining power, presented in the form 
of losses AP, is a loss of power. The greatest share in the 
losses have the AP,,, losses connected with the SP overflow 
valve. The pressure drop in the directional proportional 
valve is integrally connected with the system operation, as 
it decides the flow intensity and hydraulic cylinder speed. 
As the cylinder speed is assumed equal to v,, = 0.20 m/s and 
the load equal to Fy = 25 KN, the AP.,, losses are the second 
greatest losses in the system. Noteworthy is here much greater 
field Pw, of the cylinder useful power compared with that 
field presented in Fig. 4. This is also connected with higher 
cylinder speed. 

Also noteworthy in the presented figures are losses 
connected with the pump. Although they make a small part 
of the total losses, they are next in size to the structural losses 
AP, in the system. Mechanical losses AP,,, are in the first 
place. The field of pump volumetric losses AP,, is smaller and 
the field of pump pressure losses is smaller than that of the 
volumetric losses. 

In figures 6 and 7 the field of cylinder useful power P,,, is 
a smaller part of the field P,, of power absorbed by the pump. 
This is connected with the hydraulic cylinder small load and 
small speed. As the power of mechanical losses AP,,,, in the 
cylinder is a function of the loading force Fy, the field of power 
of those losses is also small compared with other losses in the 
system. 

Comparing the energy balance of a p = var system presented 
in Fig. 7 with the balance of a p = const system shown in Fig. 
6 with the same cylinder speed and load coefficients, it can 
be noted that the structural losses AP,, in the proportional 
directional valve and in the controlled overflow valve as well as 
the pump volumetric losses AP,,, in a p = var system are smaller 
due to the smaller pump discharge pressure pp). 

Also smaller is power AP,,, of the pump mechanical 
losses. 

Figures 8 and 9 present the energy balance of an investigated 
constant pressure (p = const) and variable pressure (p = var) 
system with the hydraulic cylinder Fu = 25 kN load and 
Vu = 0.35 m/s speed. It can be seen that structural losses are 
much smaller than in the case of limited or partial use of the 
cylinder speed and load, which is presented in figures 4 to 7. 
The pressure losses AP, in the conduits increased due to higher 
cylinder speed. 

Comparing the energy balance of the p = var system shown 
in Fig. 9 with the balance of the p = const system in Fig. 8, 
with the same cylinder speed and load coefficients, it can be 
noted that the balances are similar in the sense of areas of the 
corresponding fields. 

Structural losses are much smaller than in the case of limited 
or partial use of the cylinder speed and load. Mechanical losses 
AP vim in the cylinder and pressure losses AP, in the conduits 
increased due to higher cylinder speed and load. 


Pm 


SUMMARY AND CONCLUSIONS 


The hydrostatic drive energy efficiency is a product of the 
efficiencies of drive system components. Efficiency of the 
component elements is, in turn, a product of the mechanical, 
pressure and volumetric efficiency of those elements. While 
determining those efficiencies, the power of losses in the 
elements is not taken into account, but only torque or force 
of the mechanical losses, pressure losses in the conduits or 
intensity of the volumetric losses. It is acceptable in the case of 
rotational pumps and motors and double-piston rod cylinders. 
However, it is not sufficient in considering the energy efficiency 
of commonly used linear motors, single-piston rod cylinders 
and systems with those machines. Therefore, not always 
justified simplifying assumptions are applied [3]. 

Full picture of the energy losses in a hydrostatic drive 
system is a picture of the power of energy losses in the system 
elements. The system feeding pump shaft power is equal to 
the sum of the hydraulic motor shaft or piston rod power plus 
powers of losses occurring in the energy stream flowing through 
the component elements. Power delivered to the system on the 
pump shaft is also influenced by the interrelation between the 
pump driving motor speed np = nyn and the pump shaft torque 
Mp» Powers of the energy losses in the system elements and also 
powers developed by the elements must be precisely defined. 
The picture of energy losses requires the range to be determined 
of the hydraulic motor useful power P,,,, determined in turn by 
the range of torque M,, and angular speed ,, of a rotational 
motor shaft, or force Fy and linear speed v,, of a linear motor. 
The picture of energy losses in a hydrostatic drive system 
should be built from the hydraulic rotational motor shaft or 
linear motor piston rod towards the system feeding pump 
shaft [4, 5]. 


CONCLUSIONS 


1. The impact is presented of the Pu, power on the Pp, power 
in the analysed systems and also the impact of power AP 
of losses in the elements on the P,, power. The hydraulic 
cylinder instantaneous useful power P,,,, defined as 
a product of force F,, and cylinder piston rod speed Vy, is 
independent of all the losses. The following powers should 
be added to the useful power P u.: power AP um of mechanical 
losses in the cylinder, power AP, of losses in conduits, 
power AP „ of structural volumetric losses and power AP, 
of structural pressure losses connected with the throttling 
control, as well as powers of the pump losses: pressure 
losses AP,,, volumetric losses AP,, and mechanical losses 
AP» Lhe sum of the useful power P mw and the AP power of 
all the system losses determines the instantaneous value of 
the P,, power that the pump requires of its driving motor. 

2. The analysis of the power of losses has been performed; 
for instance, in the p = const system (Fig. 4), power APum 
of mechanical losses in the cylinder is smaller than power 
AP,,, of mechanical losses in the pump. In the p = var system 
(Fig. 5), powers of losses AP ym and AP,,,, are smaller; also 
here the AP um is smaller than the AP,,,, power. 

3. It has been shown that with cylinder unchanged load F,, 
and its increased speed v,,, the power of structural pressure 
losses AP,,, increases, because the intensity of flow through 
the proportional directional valve is increased (Figures 6 to 
9). 

4. The change of structure from p = const to p = var, with the 
same system useful power P,,,, brings in effect a serious 
decrease of the power AP,, of structural losses (Figures 8 
and 9). Simultaneously, at the same cylinder speed vy, in 
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the p = var structure, the power AP,, of volumetric losses 
in the pump and power AP,,, of mechanical losses in the 
pump decrease, but power AP,, of pressure losses in the 
pump slightly increases. 

Comparing the energy balance of the p = var system with 
the p = const system balance, with the same greatest values 
of the hydraulic cylinder speed coefficient Oy and load 
coefficient My it can be noticed that those energy balances 
are close to each other in the sense of size of the field of 
power of losses (Fig. 8 and 9). 

Considerable decrease of the power of pump structural 
losses AP,, and power of volumetric losses AP,, in the p = 
var system in relation to the p = const system can be seen 
with smaller hydraulic cylinder loads. This is connected 
with lower pp, pressure in the p = var system, because 
a system with the overflow valve controlled by the current 
pressure p, of outflow from the directional valve to the 
cylinder inlet chamber allows to adjust the pp, pressure in 
the pump discharge conduit to the current cylinder load, i.e. 
the p, pressure, in such a way that pressure loss App, in 
the directional valve slot fp,, and pressure loss App; in the 
directional valve slot fpg, are limited. The App, pressure 
drop in the directional valve is smaller than the Ap,,, drop in 
the p = const system. In connection with lower pp, pressure, 
the power Pp, absorbed by the pump is significantly lower 
(Fig. 4 and 5). 

With small hydraulic cylinder speed v,, and small load F,, 
it can be noticed, that in spite of using a constant capacity 
pump, the power of structural pressure losses AP „and also 
power of structural volumetric losses AP „is considerably 
smaller in the p = var system compared with the p = const 
system. Although the intensity of flow Q, through the 
overflow valve to the tank with the same cylinder speed 
Vy iS in the compared systems practically the same, the 
product of smaller pump discharge pressure pp, and the 
flow intensity Q, (the Q,-Q,, difference) results in a smaller 
value of the power AP,, of structural volumetric losses in 
the p = var system (Fig. 4 and 5). 

Power AP, of losses in the conduits is great at a great 
value of flow intensity Q,, towards the hydraulic cylinder 
(Fig. 8 and 9) compared with the power AP,, of pressure 
losses in the pump. In the investigated p = const and 
p = var systems, in the individual conduit sections many 
elements are installed (filters, conduit connections, cut- 
off valves in fully opened position, temperature sensors), 
which change the direction or speed of flow in the conduits. 
Therefore, even greater pump capacity Q, causes smaller 
flow resistance in the pump channels (which in the drawing 
corresponds to smaller width of the field of power AP,, of 
losses) than the flow resistance in the conduits (which in 
the drawing corresponds to the width of the field of power 
AP. of losses). 

Power AP\,,, of mechanical losses in the hydraulic cylinder 
changes with the change of pressure in its chambers. In 
a variable pressure system power AP ym of mechanical losses 
is smaller than in the p = const system, as the Fum force of 
friction losses in the cylinder is smaller (Fig. 4 and 5). 
With increasing the hydraulic cylinder speed v,, and load 
Fy to the maximum values Vmmax and Fuma respectively, the 
power AP of structural volumetric losses and power AP y 
of structural pressure losses, connected with the throttling 
control assembly, is minimized. 

When the hydraulic cylinder does not displace (when its 
speed is equal to zero - v,, = 0) and the pump operates, the 
cylinder useful power P,,, is equal to zero; the following 
powers of losses occur in the system: power AP,,, of 


stv 
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structural volumetric losses and powers of losses in the 
pump: pressure losses AP,,, volumetric losses AP,, and 
mechanical losses AP,,,. The sum of those powers is equal 
to the power P,, absorbed by the pump. With smaller loads 
Fu of the stopped cylinder, power Pp, absorbed by the 
pump in the p = var system is smaller than power P,, in the 
p = const system 

The presented energy balances indicate that power Pp. 
absorbed by the pump (with the same cylinder useful 
power Pw) is different in two investigated systems. Pump 
operating in the p = const system with constant capacity 
and nominal pressure p, absorbs all the time much greater 
power P,, than in the p = var system. Apart from the wasted 
power, elements of the constant pressure system get worn 
out faster and hydraulic oil as the working medium is also 
used faster. 

Presentation of the energy balance of investigated systems 
by means of graphical interpretation as fields of power 
of energy losses shows in a legible and simple way the 
powers in the system and powers of losses in the elements, 
dependent on the hydraulic cylinder speed v,, and load 
Fw on the speed v,, throttling control structure and on the 
quality of the system component elements. This is the first 
example of the changes of fields of power AP of losses in 
the elements and of power Pp. absorbed by the pump, as 
a function of the hydraulic cylinder useful power Pwu: 


. The variable pressure system has clearly smaller, compared 


with the constant pressure system, power of structural 
pressure losses in the throttling control assembly during 
a decreasing external load on the hydraulic linear motor. 
Also power of structural volumetric losses in the controlled 
overflow valve decreases, although the intensity of flow 
of the volumetric losses in that valve slightly increases 
compared with the constant pressure system due to greater 
pump capacity. In the pump operating in a variable pressure 
system, slight increase occurs of the power of pressure 
losses, decrease of the power of volumetric losses and also 
decrease of the power of mechanical losses. In effect, during 
a small load on the hydraulic linear motor, a decrease is 
noticed of the power absorbed by the pump from its driving 
electric motor, which, with unchanged useful power of the 
hydraulic linear motor, significantly increases the energy 
efficiency of the whole system compared with the constant 
pressure system efficiency. 


. The two compared systems may achieve, during operating 


with maximum load and maximum speed, the same 
maximum overall efficiency. The variable pressure system 
becomes then a constant pressure system and the working 
conditions of both systems are the same. 
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ABSTRACT 


Water-lubricated bearings have been applied to support ship screw shafts for over a hundered years. 
Development of plastic materials has caused that novel sliding polymers appeared available on the market 
and being approved by classification institutions, possible to be applied in shipbuilding. However in the 
subject — matter literature there is no reference to application of bearings fitted with metal or ceramic 
bushes in shipbuilding. Nevertheless they have important merits such as low resistance to motion, long 
service life or stability of form. But some doubts are evoked by their large value of shape elasticity module 
which can lead to local stress concentration at the bush edges in the case of non-axiality of bush and shaft 
axes. Hence this work has been aimed at the testing of a bearing made of bronze-graphite silter. In the 
work is presented a comparison of measurement results of resistance to motion , hydrodynamic pressure 
distributions in lubricating film as well as shaft axis trajectories of the bearing, with those made for a typical 
elastic polymer bearing. The measurement results have showed high quality of the tested material. In the 
opinion of these authors it could be applicable to bearing ship propeller shafts. 


Keywords: bearing; bearing of ship propeller shafts; tribology; water-lubricated bearings 


INTRODUCTION 


As a result of becoming screw propulsion on ships 
widespread it was necessary to elaborate a reliable and durable 
bearing for ship propulsion shafts. Most of operational problems 
have been generated by the last section of the shaft on which 
screw propeller is installed — the so called screw propeller shaft. 
It results a.o. from a complex loading system in consequence of 
mass forces and harsh working conditions, especially exposure 
to salt sea water. From the very beginning of the application of 
mechanical propulsion to ships till now a slide bearing has been 
used in the classical propulsion system of floating units. In the 
past attempts have been made to use rolling bearings, however 
such solution can be met now rather rarely as it makes repair 
technology complicated, is less durable as well as cannot be 
used for large diameter shafts [1, 2]. 

Since many years ago engineers have faced troubles 
associated with propeller shaft bearing and in the past many 
failures of shafts, bearings and sealings, have happened. With 
the famous ship Great Eastern is connected an interesting 
story dealing with her propeller shaft bearing. This unusual 
ship propelled by screw propellers and side paddle wheels, 
of more than 200 m in length, built in 1859, was then several 
times greater than any other ship. After two cross - Atlantic 
trips it was revealed that her stern slide bearing fitted with 
white metal liner was completely worn out. Measurements 
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of position of the shaft (shaft settling) showed that its wear 
reached the value of 63 mm. According to suggestion of 
John Penn, a British engineer, the white metal bearing was 
replaced by a new one made of guaiac wood (lignum vitae) 
[3]. After completing such modernization the problem of 
stern bearing did not recurred during the entire service period 
of the ship. 

Already at the beginning of 20% century such wooden 
bearings have been in common use. However they have had 
some limitations. Among the most important can be numbered 
the sensitivity to lubrication with a contaminated lubricant e.g. 
in the conditions of sailing in shallow or inland waters. 

Dynamical development of novel tribological materials 
including plastics has resulted in their successful application 
as a sliding material to water-lubricated bearings. Today there 
are a few worldwide companies which produce bearings with 
a polymer sliding bush approved by ship classification societies 
for application to ships.The tests performed in the past have 
demonstrated their high quality, hence they are widely applied 
today on both small and large ships. 

In spite of the common application of polymer bearings, 
new research projects aimed at elaboration of better and better 
sliding materials, are still under way. Bronze-graphite sinters 
of all kinds can serve as an example of a material alternative to 
plastics. They have many merits, e.g. a high stability of form, 
higher resistance to elevated temperature and seizing-up. 


ORIGIN 


For recent years many types of water-lubricated slide 
bearings have been tested at Gdansk University of Technology 
[4+8]. Among other things, have been performed tests of 
a water-lubricated bearing fitted with the sliding bush made 
of BA1032 (CuAl10Fe3Mnz2) Al-Fe-Mn bronze and journal of 
X10CrNi18-8 (AISI 301) stainless steel. The tests showed that 
such association cannot be lubricated with water, despite its 
perfect work in the conditions of lubrication with oil or plastic 
lubricant. The tested bearing has been seized up shortly after 
its starting (Fig. 1). 


Vienm im HPA) epee 


Fig. 1. The damaged shaft made of X10CrNi18-8 stainless steel-fragments 
of the broze bush can be observed on its surface 


However in recent years novel, special sliding materials 
made on the basis of bronzes, iron or nickel, capable of 
operating in water, have appeared. One of such materials is 
bronze —graphite sinter. Its very good sliding qualities in the 
case of lubrication with water, advertised by its producer, 
result from its composite structure. Into bronze-tin matrix 
a lubricating medium in the form of graphite was inserted. 
Due to this, even after breaking lubricating film graphite 
particles effectively lower friction forces and prevent seizing. 
The material can be applied even in the temperature close 
to 800°C and is of a high resistance against corrosion. 
As novel sintered materials are successfully applied to water- 
lubricated bearings of vertical water turbines it was decided 
to test a material of the kind in the operational conditions 
of ship stern tube bearing. It should be stressed that screw 
propeller shaft which operates in horizontal position, produces 
much more harsh operational conditions for such bearing in 
comparison with those for a vertical shaft. It results from greater 
radial loads due to statical and dynamical mass forces occurring 
in the shaft - screw propeller system. 

It was decided that such metal bearing will be compared 
with a typical, water-lubricated polymer bearing approved by 
classification institution. Such comparison will make it possible 
to determine drawbacks and merits of both the solutions. 

For the tests was selected a bronze-graphite silter offered 
in the form of prefabricated bush elements pressed under high 
pressure and temperature. The bronze bush was pressed into 
a steel sleeve and then the internal hole was machined to its 
final dimensions. For the reason ofa relatively high hardness of 
the implemented bronze (in comparison with that of polymer) 
the bush of the kind has been intended for interaction with the 
journal of the hardness of 35HRC at least. Therefore on the shaft 
made of X10CrNi18-8 steel a hardened-tempered sleeve made 
of X30Cr13 (AISI 420) steel was bonded by using thermo- 


compression method. The sleeve was next grinded and polished 
to ensure appropriate coaxiality with respect to journals of 
bearings which supported the shaft during the tests. 


a 


Fig. 2. The tested bearing fitted with metal bush: a) the prefabricated 
sliding bush element made of siltered bronze, b) the complete bearing 
before its mounting at the test stand, in which the sliding bush 
and outer steel sleeve with circumferentially placed sockets 
for pressure transducers, are seen. 


PLAN OF THE TESTS 


Data of the tested bearings are presented in the tab. 1. 

As resulted from measurements, the polymer bearing had 
a non-ideal form (cylindricity error) as well as a worse state 
of sliding surface in comparison with that of the bearing made 
of silter. It results from difficulty in machining a very elastic 
polymer and its deformations resulting from absorption of 
water (soaking). 

Many calculations performed in the past, based on 
hydrodynamic lubrication theory for EHL model [5, 6, 7] have 
shown that the bearing with rigid bush operated at a larger 
film thickness in comparison with the same bush but made of 
polymers. The different thickness of lubricating film in rigid 
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Tab. 1. Comparison of the tested bearings 


Shaft diameter 
[mm]/Bush 
length [mm] 


Shear modulus 
[MPa] 


Material 


Bronze-graphite 
silter 


Uniform elastic 
polymer approved 
for marine 
applications 


100/150 


bearing and elastic one results from elastic bush deformation 
under pressure, which changes lubricating gap form [9]. 


= 
o 


--m- -polymer bush E=600MPa 
~—®-~ polymer-composite E=4000MPa 
—4— brass E=10000MPa 


average pressure in hydrodynamic film [MPa] 
o o ° o 
N > o œ 


4 
© 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
diameter clearance [mm] 


Fig. 3. Influence of shear modulus of bush material on load-carrying 
capacity of bearing, in function of bearing clearance value; journal 
diameter of 100 mm, shaft rotational speed of 11 rev/s 


One of the main aims of experimental tests was to compare 
theoretical results with experimental ones as well as to determine 
drawbacks and merits of the bearings with metal bush. 


EXPERIMENTAL TESTS 


For the experimental tests was used the test stand already 
described in the subject-matter literature [6,10]. It has wide 
research possibilities. During the tests, resistance to motion, 
pressure in lubricating film as well as shaft axis trajectory was 
measured. 

The tests were performed for working parameters typical for 
marine application of such bearings. In order to gain as much 
information about tribological qualities of both the bearings 
as possible,was planned a set of tests in different conditions 
realized by changing both load and rotational speed. The load 
was applied statically by hanging a weight. The tests were 
performed for three different load values giving the average 
pressures of 0.2, 0.4 and 0.6 MPa, respectively. Rotational speed 
was varied in the range from 0 tol! rev/s. The measurements 
were taken for the speed values of 1, 2, 3, 5, 7, 9 and 11 rev/s. 
In the case of determination of resistance to motion also statical 
friction coefficient was measured. 

During the tests, water of 15°C temperature was pressed 
through the bearing. Its flow rate amounting at least to 6 l/min 
was sufficient to ensure proper lubrication of the bearing. All the 
measurements was performed for the initially run-in bearings 
(after 24h working period in the range of low rotational speeds, 
at mixed friction work). 
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Bearing 
clearance [mm] 


Sliding surface 
roughness 
R, [mm] 


Shaft journal 
roughness 
R, [mm] 


0.17-0.19 
(shape errors) 


RESULTS OF THE TESTS 
Tests of starting resistance 


A very important feature of water-lubricated slide bearings 
is their capability of starting under load. This is especially 
important in the case of horizontal ship propeller shafts which 
exert their weight onto bearings regardless of whether the 
shaft is in idle or motion state. In slide bearings, in contrast to 
rolling ones, the problem of appearance of very large values 
of resistance during start-ups under load, often happens. It 
results from a large, statical friction coefficient which occurs in 
conditions of lacking hydrodynamical lubrication. In the case 
of polymer bearings the effect can be heightened by stick-slip 
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Fig. 4. Measured friction coefficient in: a) bearing with metal bush, 
b) bearing with polymer bush 


phenomena. Usually, after reaching even a low rotational speed 
(of abt. a few rpm), a partial raise of hydrodynamical pressure 
occurs and in consequence friction coefficient drops suddenly. 
In Fig. 4, are presented the diagrams of measured friction 
coefficients in metal bearing and polymer one in function of 
rotational speed and average pressure values. 

The recorded friction coefficient values indicate that in the 
bearing of metal bush a larger value of friction coefficient equal 
to abt. 0.5, occurs. In the same conditions in the polymer bearing 
the friction coefficient amounts to abt. 0.25. It means that the 
bearing with metal bush requires at least two times greater 
starting moment in comparison with the polymer bearing. It 
is also seen that in the metal bearing dropping resistance to 
motion is faster because after reaching the rotational speed 
of 0.5 rev/s the friction coefficient drops to the five times 


lower value equal to abt. 0.1. It is undoubtedly connected 
with forming hydrodynamic lubricating wedge in the rigid 
bearing faster than in the elastic (polymer) bearing. It should be 
mentioned that in the bearing with metal bush lower resistance 
to motion occurs in steady working conditions (regardless of 
pressure values). 


Tests of journal displacement against bush 


The trajectory diagrams are presented on the background 
of experimentally measured clearance circles (Fig. 5 + 10). It 
makes a better interpretation of the results, possible. Moreover, 
shape of the clearance circles reflects shape of the bush, hence 
it makes it possible to assess bush shape errors. The shape error 
which appears in the bearing of polymer bush is not surprising 
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Fig. 5. Shaft orbits measured on both sides of the bearing with metal bush, under load of 0.2 MPa; 
a) left side of the bush, b) right side of the bearing; orbits for three successive shaft revolutions 
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Fig. 6. Shaft orbits measured on both sides of the bearing with metal bush, under load of 0.4 MPa; 
a) left side of the bush, b) right side of the bearing; orbits for three successive shaft revolutions 
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Fig. 7. Shaft orbits measured on both sides of the bearing with metal bush, under load of 0.6 MPa; 
a) left side of the bush, b) right side of the bearing; orbits for three succesive shaft revolutions 
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Fig. 8. Shaft orbits measured on both sides of the bearing with polymer bush, under load of 0.2 MPa; 
a) left side of the bush, b) right side of the bearing; orbits for three succesive shaft revolutions 
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Fig. 9. Shaft orbits measured on both sides of the bearing with polymer bush, under load of 0.4 MPa; 
a) left side of the bush, b) right side of the bearing; orbits for three successive shaft revolutions 
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Fig. 10. Shaft orbits measured on both sides of the bearing with polymer bush, under load of 0.6 MPa; 
a) left side of the bush, b) right side of the bearing; orbits for three successive shaft revolutions 


as it results from machining the elastic water-absorbing 
polymer. However, is puzzling the shape error of the bearing 
with metal bush machined to its final form after installation 
into the metal sleeve,. The tests of the bearing made of BA1032 
(CuA110Fe3Mn2) bronze, whose sliding surface was grinded, 
carried out in the past (as already mentioned), showed that it 
is possible to obtain its shape almost ideal (Fig. 11). In the 
discussed case the bush length was as large as 300 mm and in 
spite of that the high machining precision (H7 tolerance class) 
was achieved. 

As results from the shaft axis trajectory measurements, the 
bearing operated at a significant non-coaxiality of shaft against 
bush (Fig. 5 + 10). In spite of that it operated in the range of 
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fluid friction. In the opinion of these authors it testifies to their 
high potential and big margin of load-carrying capacity. It 
should be added that the revealed non-coaxiality does not result 
from inappropriate assembling work but from free-support of 
the bearing, hence even a small asymmetry in hydrodynamic 
pressure distribution produces natural twisting of the bush. 

Repeatability of trajectories draws attention. In the diagrams, 
orbits of three successive revolutions of shaft are presented. In 
most of the cases the oval orbits smaller than 10um in size cover 
practically each other. Due to this the bearing practically does 
not produce vibrations, except of the situations of its unstable 
work (Fig. 8). Such phenomenon appeared only in the bearing 
with polymer bush under load of a small value. 
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Fig. 11. Experimentally measured clearance circles of BA1032 bronze bearing whose sliding surface was grinded; 
a) left side of the bush, b) right side of the bearing 


Tests of hydrodynamic pressure distribution 


Below in Fig. 12 + 14, are presented selected results of 
experimental tests of hydrodynamic pressure distribution in 
the metal (rigid) bearing and the polymer (elastic) one. By 
analyzing the obtained results the significant effect of bush 
rigidity on pressure distribution form can be observed. In the 
rigid (metal) bearing much greater maximum pressure which 
varies in function of values of average pressure and speed, 
occurs. The averaged values of the results achieved during the 
tests are given in Tab. 2 below. In the case of the rigid bush the 


maximum pressure is up to four times greater than the average 
nominal one, whereas for the flexible bush the maximum 
pressure does not exceed twofold value of the average pressure. 
It results from local change of lubricating gap form in the 
flexible bush, that leads to dropping maximum pressure at 
simultaneous increasing lubricating wedge lap angle. 

During operation of the bearing in normal conditions 
pressure pulsations are rather small (Fig. 15a). As expected, 
the above described unstable work of the bearing (Fig. 8) 
generates significant pressure pulsations in the lubricating 
film (Fig. 15b). 


Tab.2. Average values of pressure in lubricating film 


Shaft rotational 
speed [rev/s] 


Nominal (average) pressure 
resulting from load [MPa] 


Maximum pressure in 
lubricating film [MPa] 


Rigid (metal) bush 


Flexible (polymer) bush 
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Fig. 12. Measured distributions of pressure in lubricating film in lower half of the bush under load of 0.2 MPa; 
a) bearing with metal bush, b) bearing with polymer bush 
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Fig. 13. Measured distributions of pressure in lubricating film in lower half of the bush under load of 0.4 MPa; 
a) bearing with metal bush, b) bearing with polymer bush 
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Fig. 14. Measured distributions of pressure in lubricating film in lower half of the bush under load of 0.6 MPa; 
a) bearing with metal bush, b) bearing with polymer bush 
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CONCLUSIONS 


e Analyzing the measurement results one can observe distinct 
difference between behaviour of the rigid (metal) bearing 
and flexible (polymer) one. Both the bearings generate 
hydrodynamic pressure distribution of different form. 
In the rigid bearing the maximum pressure is up to four 
times greater than the average nominal pressure, and in the 
flexible bearing — two times greater at the very most. In the 
metal bearing the starting moment over two times greater 
than that in the polymer bearing, occurs. Statical friction 
coefficients amount to 0.5 and 0.25, respectively. Worth 
adding that in the polymer bearing stick-slip phenomena 
were observed, that could detrimentally affect its friction 
coefficient value. 

e Hydrodynamic pressure very fast develops and the friction 
coefficient value recorded during the tests for operation in 
the fluid friction range, amounted to 0.01, on average. It 
should be remembered however that due to limitations of the 
test stand the value is burdened with an error and evidently 
greatly overstated. 

e The water-lubricated bearing with steel bush is of a greater 
load-carrying capacity in comparison with a typical polymer 
bearing. It results from lack of deformations of the rigid 
bush under hydrodynamic pressure as well as a high 
smoothness of its surface. This way the experimental tests 
confirmed the calculation results (Fig. 3). 

e The bearings with steel bush show a lower sensitivity to 
contaminations which do not stick in sliding surface. Due 
to this fact wear occurs mostly on the side of bush, which 
is favourable from the point of view of repair cost. 

e Apparently the necessity of application of hardened stainless 
steel journals was an obvious disadvantage of the metal 
bearings, which greatly rises cost of realization of such 
solution. 


RECAPITULATION 


The performed tests showed that in the assumed working 
conditions the water-lubricated bearings with metal bush are 
of a sufficient load-carrying capacity and can be applied also 
in the case of ship propeller shafts. Approval of the applied 
material by ship classification institutions remains to be 
a formal problem only. 

The producer of the material indicates that one of the 
cardinal merits of such designs is their much greater resistance 
to wear in the conditions of lubrication with contaminated 
water since in the metal bearings, in contrast to polymer ones, 
the contaminations do not stick into bush sliding surface and, 
consequently, do not scratch the journal. As the referred opinion 
evokes some doubts, these authors have attempted to prepare 
tests of bearings of the kind in the conditions of lubrication 
with contaminated water. However the bearing with metal 
bush has an important disadvantage. It should be remembered 
that in the case of ship main propulsion shafts a problem is 


their appropriate assembling, and as a rule we then deal with 
a smaller or greater non-coaxiality between shaft and bush. 
It may also occur occasionally as a result of deformation of 
a flexible ship hull at its non-uniform weight distribution of 
cargo. The effect may also be strengthened due to hull response 
to waves. Therefore in some cases the so called edging, i.e. 
concentration of pressure between shaft journal and bush 
edge, may happen. The phenomenon can generate sudden 
wear process of the sliding pair. It will not occur in the case of 
applying elastic polymer because the bush will be deformed 
elastically, that leads to a greater contact zone. 
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Multi-stage magnetic-fluid seals for operating 
in water — life test procedure, test stand and 
research results 
Part I 
Life test procedure, test stand and instrumentation 
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ABSTRACT 


The edge criteria for magnetic-fluid seals durability research are critical for the seal 5 
technological usefulness. The magnetic fluid seal technology in water is much more 
difficult technological problem in comparison to gas or vacuum environment. There is 
always problem with immediate contact between magnetic fluid and the environmental 
uninspected liquids. In the paper are presented ways of experiments definitions carried out 
with multi-stage rotating magnetic fluid seal operating in permanent contact with utility 
water. The special test procedures have been elaborated and practically used, with the 


aim to define main characteristics of magnetic fluid seals. These characteristics were: critical pressure, 

critical motion velocity and working life of the seal. Some grades of commercial, silicon based magnetic 

fluids of hydrophobic properties were tested in sealing system. The test stand and its features as well as 

results of tests are presented in the form of tables, cross section draws and diagrams. In conclusion one 

can see that multi-stage magnetic-fluid seals could be efficiently used in water for rotating shaft seals in 
a limited range of motion velocity and work lasting. 


Keywords: shaft sealing; utility water; seal durability; critical pressure; critical motion velocity; magnetic fluid 


INTRODUCTION 


The investigations described in this paper were performed 
at the test stand designed and built in accordance with the 
below presented assumptions, especially intended for testing 
drives used in ocean engineering. The research was finacially 
supported in the frame of research projects of Ministry of 
Science and Higher Education and NCBR. Many elements 
of the research methods, measurement data recording and 
maintenance techniques are entirely original. Therefore it 
is necessary to describe the applied research process and 
procedures. Results of testing the single-stage magnetic-fluid 
(MF) seals have been already published earlier and those 
dealing with multi-stage ones are given in the second part of 
this elaboration. The problem in question is below presented 
on the background of the state of engineering with reference to 
the most important details given already in the paper [18]. 

This paper relates results of a continuation of the tests of 
MF seals in liquid environment, carried out in 2011, this turn 
with the use of multi-stage sleeves. Worth reminding that 
the magnetic fluid seals operating in liquid environment are 
exposed to loss of tightness due to action of factors which occur 
during direct inter-phase contact of two environments being 
in the same physical state — magnetic fluid and sealed liquid. 


62 POLISH MARITIME RESEARCH, No 4/2012 


It is hard to compare the problem with the situation where gas 
environments have to be separated. The basic condition to be 
fulfilled in designing the seal is to select a magnetic fluid (or 
environment liquid) of such physical chemical properties that 
the liquids being in contact to each other do not mix together 
and on their separation surface a high inter-phase surface 
tension occurs. In statical conditions (i.e. when seal elements do 
not displace against each other), if only the system is properly 
designed, satisfaction of the above mentioned conditions 
usually ensures correct work of the seal. Lack of susceptibility 
to mixing and high inter-phase surface tension are as arule 
insufficient if a MF seal has to be applied in conditions of seal 
motion or sealed liquid flow in vicinity of the seal. Relative 
motion of the two liquids generates a surface instability on the 
boundary of the phases [2, 3], whose intensity greatly increases 
above certain limiting difference in velocity, that results in 
washing away the magnetic fluid and loss of tightness of the 
system along with time. Life time of such seal drops below its 
profitability. 

In the published elaborations [4, 5, 6] are described attempts 
at determining possible conditions of operation of MF seals in 
water, depending on pressure and motion speed. In the work 
[7] are presented investigations on application of MF seal to 
sealing lubrication oil, and in [8] is given an attempt at design 


optimizing the seal intended for the operating also in oil. In the 
paper [9] its authors present a concept of the „magnetic surface 
tension” and its impact on MF seal operation in water. 

Some years earlier in AGH Laboratory, the tests of MFseals 
operating in utility water environment were carried out [12, 
13, 14] in order to assess their applicability in shipbuilding. 
Results of the tests conducted in the years 2009+2010 on the 
test stand adjusted to operation of the seals in liquids [15] 
showed that it is possible to apply some commercial magnetic 
fluids to MF seals operating in direct contact with water in 
working conditions which occur in many devices met not 
only in ocean engineering. In the tests of a utilitary character 
single-stage MF seals — models of real seals with magnetic 
fluid - were used; their working conditions were limited only 
to the assigned technical task, i.e. assessment if performing 
the tests is purposeful. 

In this work are described experimental tests of a magnetic 
fluid seal of three sealing stages, carried out in its direct contact 
with water and in ship propulsion conditions as well. 

Alike as at the tests of the pack with single-stage sleeve, 
the tests were performed in variable magnetic field conditions 
(different number of permanent magnets in magnetic system of 
the seal), for different set-up values of working pressure and the 
same magnetic fluid dose delivered to the sealing stages. 

In the presented tests a sleeve of three unsymmetrically 
arranged sealing stages was used. The application of the 
unsymmetrical (one-sided) arrangement of sealing stages had 
the practical aim similar to that as in the tests of single-stage 
seal [18]. In such arrangement the closing of magnetic field is 
realized by the cylindrical part of the sleeve fitted with a special 
longitudinal groove, that makes it possible to drain away any 
leakage and diagnose continuously working state of the seal. 

For realisation of the tests was built a test stand, called 
shortly MFSL-135, composed of a special, purpose-oriented 
testing head, driving system containing an electric motor 
and its operation controller, as well as measuring system 
composed of gauges and transducers for measuring pressure, 
torque of driving motor and detecting leakage from the tested 
seal. The test stand is fitted with multi-channel measuring 
aparature for measurement, data acquisition, processing and 
recording. 

In the research concept it was assumed to perform short and 
long - term measurment tests which would make it possible 
to determine limiting operation features of the MF seals in 
water, and consequently to predict service life of seals of the 
type operating in working systems of ocean engineering units. 
To check influence of both external and structural factors on 
the seal’s operation, it was assumed to perform the tests under 
different working pressure values of water as well as with 
switched-on magnetic systems adjusted to generating magnetic 
field of different intensity within the sealing system. 

The design solution of the testing head and measuring 
system enables also to perform tests of other kinds, e.g. tests of 
hybrid seals composed of a magnetic-fluid seal and pneumatic 
shield which prevents the first from direct contact with water. 
The tests of the MF seal with pneumatic shield performed in 
the frame of this project were presented in the form of a lecture 
given during a technical scientific conference [19]. 

In this paper are presented: the test stand, research concept 
and assumptions, as well as the criteria necessary for assessing 
correctness of measurement tests. Also, a programme of the 
tests, their scope and instrumentation as well as magnetic fluids 
used for the tests, are described. Moreover, testing procedures, 
technique of preparation and performance of the tests are 
highlighted in this part, and the research results are presented 
in the second part of the article. 


RESEARCH CONCEPT 
Research hypotheses 


Effectiveness of MF seal in the case of direct contact 
between ferromagnetic fluid and water on the inter-phase 
surface depends on many external factors as well as physical 
chemical features of the ferromagnetic fluid used in the seal. 

On the basis of the observations and research experiments as 
well as determination of magnetic fluid features, made during 
the tests of single-stage MF seals [18], was elaborated the 
research concept adjusted for conducting the tests of multi-stage 
MF seals operating in water, and assessing their results. 

In the tests in question were used magnetic fluids of 
physical chemical features determined by their producer, and 
standardized testing procedures. Elaboration of appropriate 
procedures and obeying them was a necessary condition for 
correct assessing influence of external factors on effectiveness 
of seal operation. 

The following basic research hypotheses concerning 
influence of external factors were advanced: 
¢ Hypothesis I 

- Life of MF seal depends on relative speed of its motion 

in utilitary water environment. 
°- Hypothesis II 
- In the critical pressure range specific for a given MF 
seal its life depends on relative value of set-up working 
pressure in relation to critical pressure. 


Relative pressure 


Loss of tightness of multi-stage MF seal is a much more 
complex process than that in the case of single-stage seal. It 
consists in penetration, resulting from exceedance of critical 
pressure in successive sealing stages; total penetration pressure 
of multi-stage seal is asum of pressure increments causing 
penetration of successive stages. Worth mentioning that in most 
cases critical pressure of multi-stage seal is lower than simple 
multiplicity of critical pressure of single-stage seal of the same 
geometry and magnetic field parameters [20, 21]. However for 
purposes of the realized tests was made a simplified assumption 
that penetration pressure values of successive stages are equal 
to each other. 

As absolute value of working medium pressure greatly 
affects seal life [18] the original concept of relative pressure 
was used in the research on magnetic fluids of different critical 
pressure values. With a view of complexity of the problem the 
concept in question should be highlighted. 

In this research, for each of the tested MF fluid, critical 
pressure tests were performed with the use ofa three-stage seal. 
The critical pressure of one stage of the seal can be defined 
as follows: i 


Pki = 3 Piet (1) 


On the basis of the assumed critical pressure of one stage 
of the seal the relative pressure was defined as the ratio of the 
set-up working pressure of the seal and the critical pressure 
determined for one sealing stage: 


Py= t= 3 Pe D 
Pkrri Pkrin 
where: 
Pp, — relative pressure, 
p. — working pressure (set-up for test), 
Part critical pressure of one sealing stage, 
Pam total critical pressure of a three-stage seal. 
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RESEARCH ASSUMPTIONS 


The following research assumptions were made: 

e Assumption I 
In the tests magnetic fluids of hydrophobic properties, 
produced by specialized producers (FerroLabs) [17], 
available on the market, will be used. 

e Assumption IT 
The tests will be conducted with the use of three-stage 
seals of 50 mm nominal diameter, already implemented in 
the previous tests, and 0.1 mm height of the gap filled with 
magnetic fluid, that should make it possible to compare 
results of these tests and those performed earlier. 

e Assumption II 
Magnetic field which keeps magnetic fluid in the gap, will 
be generated by the set of permanent magnets in the form 
of discs circumferentially and uniformly located between 
pole shoes of magnetic system. The tests are assumed to be 
performed for different magnetic fields, that will be realized 
by changing number of magnets in the magnetic system of 
the seal. 

e Assumption IV 
The highest motion speed at which, after reaching 1 mln 
cycles (number of rotations), the tested seal does not yet 
lose its tightness, will be called the long-term operation 
speed of a seal in given working conditions. 


MEASUREMENTS AND CRITERIA FOR 
THE DETERMINING OF LIMITING 
PARAMETERS 


Measurements of limiting parameters 


The measurements were conducted on the basis of standard 

measurement tests of three types. 

¢ The test of total statical critical pressure of the seal. 
The test is realized in statical conditions (motionless shaft). 
During the test the pressure is continuously increased till 
occurrence of seal leakage symptoms. The test is aimed at 
determination of a pressure value which will serve as the 
basis for the determining of relative pressure values to be 
applied in the tests. 

¢ The test of critical speed of the seal 
The measurement of seal critical speed results in continuous 
increasing rotational speed of the shaft with a constant 
acceleration, the same for all the tests, till occurrence of seal 
leakage symptoms. For each relative pressure value used 
in the research the separate critical speed test is performed. 
The test is aimed at determination of the highest speed 
beginning from which the series of measurements of long- 
term operation of the seal is commenced. 

¢ The test of long-term operation of the seal 
The test of long-term operation of the seal is composed 
of the series of measurements conducted for a set-up 
relative pressure and constant rotational speed of the 
shaft. Successive measurements are taken at lower and 
lower rotational speed till reaching the speed at which no 
leakage of the seal is observed during the period of 1 mln 
cycles of its operation. On the basis of the obtained results 
the diagram of seal life versus its motion speed at a given 
relative pressure, is elaborated. 


Criteria and way for the determining of values of 
limiting parameters 


For all the tests lack of leakage was the criterion of correct 
operation of the seal. 
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¢ The total statical critical pressure of the seal, Pirr 
The total statical critical pressure is such pressure value at 
which the first leaking drop is observed in the drainage tank 
or a leakage is recorded by the leakage detector. 

¢ The statical critical pressure of one seal stage, Pyr 
The statical critical pressure of one seal stage was assumed 
to be equal to one third of the total statical critical pressure 
of the seal: Piri = 1/3 Pian 

e Critical speed of the seal 
The critial speed of the seal is the lowest value of seal 
motion speed at which the first drop of leakage is observed 
in the drainage tank during measurements with continuously 
increasing speed. 

¢ Long-term operation speed of the seal 
The long-term operation speed of the seal is equal to the 
highest rotational speed of the seal, at which no leakage of 
the seal is observed during the period of 1 mln cycles of its 
operation. 


SCOPE AND PROGRAMME 
OF THE TESTS 


Scope of the work to be performed during the tests in 
question comprised both short-term and long-term tests. In 
the tests conventionally called shot-term the statical critical 
pressure and critical speed of the seal was determined, whereas 
the tests performed in an extended period, conventionally 
called long-term ones, were aimed at determination of seal 
life values. 

The measure of pressure, used in the tests, was the relative 
pressure p,/p,,,, determined as the ratio of the set-up working 
pressure and the critical pressure of one seal stage, highlighted 
and calculated according to the principle given in Ch. 2.2. 

In spite of that occurrence of leakage was the main criterion 
for assessing correct operation of the seal, the following 
was also measured during the tests: friction moment, seal 
temperature, shaft rotational speed and pressure. The quantities 
were used to control run of the tests. 

The main tests were performed with the use of magnetic 
fluids of three kinds applied to the seal of three stages with 0.1 
mm magnetic fluid gap height. 

In the frame of supplementary research the tests were 
performed with one kind of magnetic fluid and different 
numbers of permanent magnets in the seal magnetic system. 


Programme of the main tests 


The following programme of the tests was prepared: 

e Tests of the total statical penetration pressure of the seal. 

e Tests of the critical speed of seal operation at the values 
of the relative pressure p,/p,,; equal to: 0.8; 1.6 and 2.4 
multiciplicity of the one-seal-stage critical pressure. 

e Tests of the seal life at the values of the relative pressure 
P/Pkn equal to: 0.8; 1.6 and 2.4 multiciplicity of the 
one-seal-stage critical pressure and different motion speeds 
beginning from that equal to or near the critical seal speed 
up to that at which the seal life reaches at least 1 mln 
cycles. 


In the main tests the three selected kinds of magnetic fluid 
were tested. The performed measurements were later repeated. 
If the result of the first test differed from the repeated one by 
more than 10 %, the next, explanatory, test was performed. 
During running the tests some of them were repeated if 
necessary. 


Extensometric beam 


—_— 


` . i 
Strain gauge 


Fig. 1. Extensometric leakage detector at the test stand, in its working state 


Programme of the supplementary tests 


The supplmentary tests were performed with the use of 
one magnetic fluid (FLA 003.45), selected for these tests on 
the basis of the results obtained from the main tests. The tests 
in question were conducted with the number of permanent 
magnets decreased to 10 pieces and the number of them 
increased to 18 pieces as compared with the permanent magnets 
of 14 in number used in the main tests. 

In the future is planned realization of the programme of 
supplementary tests in the same range as that above mentioned 
for the main tests, as well as determination of minimum 
intensity of the magnetic field related to N35 standard magnets 
(of the same volume). 


MFSL-135 TEST STAND 


The test stand symbolically marked MFSL-135, was 
designed and built for purposes of this research project. It 
serves for the testing of magnetic fluid seals operating in contact 
with water or another non-aggressive liquid, and may be also 
used for the testing in the environment of air or another non- 
aggressive gas. Dimensions of the stand’s head allow to place 
in it magnetic fluid seals of up to about 80 mm in diameter. For 
filling the gaps of sealing stages it is necessary to use only the 
small amount of magnetic fluid, i.e. of abt. 50-200 ul. 

The MFSL 135 test stand makes it possible to determine 
statical seal penetration pressure, limiting seal operation speed 
and seal life at a given rotational speed in the range up to 100 
s- and pressure up to 0.5 MPa. 


Design of MFSL-135 test stand 


The multi-functional test stand MFSL-135 intended for the 
testing of magnetic fluid seals was presented in detail in the 
previous papers published in this journal, and in the papers [1] 
and [18].The full list of its particular units and elements was 
also given there. 


Leakage detector 


During the previous tests, necessity of automatic leakage 
detection was stated in connection with troubles in unambiguous 
identification of instant of a leakage. Leakage identification 
during the tests with the use of the MFSL 135 stand was initially 
based on observation of outflow of water drops from the testing 
head outflow hole since drop of water pressure or disturbances 
of torque values during operation of the three-stage seal did not 
provide an appropriate diagnostic signal. As during the long- 


term testing it was not always possible to observe a leakage, 
an automatic leakage detecting system was designed. 

In the photo, Fig. 1, presenting the test stand the leakage 
detector located in the outflow zone of leaking drops, was 
shown. 

The device is composed of a high-sensitive extensometric 
beam fastened in a holder, with a glued strain gauge and 
a falling plate located just under outflow of drops from the 
casing of MFSL-135G head. The strain gauge was connected 
to an amplifier and one of the measuring card channels, that 
ensures recording even single drops flowing out of the testing 
head casing. 


Measuring system of MFSL 135 test stand 


In the paper [18] the schematic measurement diagram of 
the MFSL 135 test stand was shown. The list of particular 
units and elements of the stand was given in the description 
of the figure. 


Design of particular elements 
of MFSL 135 test stand 


The structural solid model of the test stand was presented 
in the previous papers published in this journal, e.g. [18]. 


MFSL-135G testing head 


In Fig. 2 is shown the cross-section of the MFSL-135G testing 
head, where its most important parts are shortly named. 


coolant inlet 


/ pole shoes 


unit of slinger discs 


pole shoes casing 


pressure measuring 
channel 


pressure chamber 
casing 


.__ testing sleeve 


pressure chamber / permanent \, coolant outlet 


magnet 


Fig. 2. Model of MFSL-135G testing head, elaborated 
by using Solid Works software 
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The shaft made of a non-magnetic steel is rotationally seated 
in the casing by means of the unit of rolling bearings: radial 
and angular one. At the end of the shaft is placed a replaceable 
testing sleeve made of a steel of high magnetic permeability. 
Axial position of the testing sleeve is fixed by means of the cap 
fitted with internal threaded hole. The cap, after being screwed 
on the threaded end of the shaft, presses down the testing sleeve 
to the unit of slinger discs, and then to the rotary shaft shoulder. 
The testing sleeve is placed in the hole of the unit of pole shoes 
made of a material of high magnetic permeability, seated in 
the casing of a non-magnetic material. Between the pole shoes 
the unit of permanent magnets is located. The axial position 
of the pole shoes is fixed by means of the pressure measuring 
disc and pressure chamber casing which is screwed down to 
the pole shoes casing by using the unit of bolts and nuts. In the 
pole shoes casing the inlet and outlet holes are made for the 
medium which cools the pole shoes during testing. 

In Fig.3 the three-stage testing sleeve used in the tests in 
question, is presented. In the part a) of the figure the cross- 
section of the sleeve solid model is shown, and in the part b) 
— the view of the testing sleeve after series of the performed 
tests. 


a) A b) 


Fig. 3. Three-stage sleeve: a) model; b) real sleeve 


In Fig. 4 is presented design structure and operation 
principle of the tested magnetic-fluid sealing system of three 
sealing stages. It shows the course of magnetic field lines, 
location of the sealing stages and cylindrical part with leakage 
draining groove, as well as location of water zone in the testing 
chamber. Direction of water flow through the leakage 

draining grove after loss of tightness by the sleeve, is 
indicated by the arrow. In the system the rotary slinger disc is 
placed between the sleeve with stages and the shaft shoulder, 
and the permanent disc - between the unit of pole shoes and 
the rear wall of the pole shoes casing. 


seal with magnetic fluid 
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water cooling 


water supply 
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~ magnetic 
a field lines 
rove 4 three-stage non-magnetic materials 


magnetic fluid sleeve 


Fig. 4. Design structure of the tested magnetic fluid sealing system with 
three sealing stages 
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The rotary disc is aimed at discharging leaking drops 
onto internal surface of the casing from which they are 
gravitationally discharged through the leakage outflow hole 
located in the lower part of the casing, down to the leakage 
detector plate to trigger recording the leakage instant by 
LabView software. 

In Fig. 5 are shown details of way of bearing the shaft of 
MFSL-135G testing head as well as of tension control of the 
angular bearing. From the side of the pole shoes the shaft is 
seated in the radial bearing whose axial position is fixed by the 
internal retaining ring. On the shaft, from the side of the drive, 
is mounted the angular bearing whose tension is controlled by 
means of the threaded bearing sleeve. 


pole shoe pole shoes casing 


. 


pole shoes casing 


pole shoe 


shaft controlled 


radial bearing pressing down 


shaft 
—~ rotary disc 
motionless 
disc , angular 
leakage outlet L discs radial bearing | bearing 


Fig. 5. Slinger discs: motionless and rotary one, 
mounted in MFSL-135G head 


GEOMETRICAL AND STRUCTURAL 
CHARACTERISTICS OF THE SEAL 


In Fig. 6 are shown basic dimensions of the sleeve with 
sealing stages, arrangement and dimensions of the stages, as 
well as dimensions of the magnetic system of the seal. 


52e 
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Fig. 6. Basic dimensions of the sleeve with sealing stages and dimensions of 
the magnetic system of the magnetic fluid seal used in the tests 


In Tab. 1 are collected the geometrical and structural 
parameters ofthe magnetic fluid seal. The parameters are divided 
into changeable and constant ones during testing cycle. 


MATERIALS AND DEVICES 
USED DURING THE TESTS 


Magnetic fluids 


The tests were conducted with the use of three magnetic 
fluids characterized in Tab. 2. In the table physical features 
of the tested fluids are compared and their typical technical 
applications presented. 


Tab. 1. Geometrical and structural parameters of the magnetic fluid seal used in MFSL-135G testing head 


PARAMETER 


Quantities constant during testing 


Shaft diameter, d, 


External diameter of pole shoes, d, zew 


Diameter of three-stage sleeve, dyn 


Number of stages of testing sleeve, z 


Nominal gap height, 6 


Diameter of permanent magnets, d,, 


Height of permanent magnets, h 


Amount of magnetic fluid delivered to each of the stages 


Number of permanent magnets fastened in the head casing 


10; 14; 18 


Working medium 


water 


Quantities changeable depending on a type of testing 


Rotational speed 


1— 100 rps 


Range of set-up values of working pressure 


0 — 0.22 MPa 


Set-up values of relative pressure, p,, = P,/Pkri 


P/P = 0.8; 1.6; 2.4. 


Cooling intensity 


Flow controlled by indications of head temperature 
transducer 


Head starting-up duration time up to instant of reaching the 


working speed 


250 s — for measurements of limiting speed 
20 s — for measurements of operation time to 
leakage 


Tab. 2. The most important parameters of the tested ferromagnetic fluids [17] 


PARAMETER 


Kind of fluid: 


FLS 040.040 


Kind of suspension vehicle: Silicone liquid 


FLA 002.25 FLA 003.45 
Siloxane Siloxane 


Saturation magnetization, A/m 40 + 50 


25 45 


Working temperature, °C -70 do 150 


-70 do 150 -70 do 150 


Critical temperature, °C 200 


250 200 


Freezing temperature, °C -100 


Lack of data Lack of data 


Plastic viscosity in 27°C 


0.3-0.8 Pa's 
temperature 


700-750 mPa:s 400-500 mPa:s 


Colour of fluid: black-brown 


black-brown black-brown 


Evaporation index, 


g- Sm2c- Lack of data 


1.0- 107 1.0 - 107 


Heat conductivity coefficient 


in 38°C temp., mWt/m-K Lack of data 


150 150 


Applications: 


Permanent magnets 


In the tests the permanent magnets in the form of discs of 
ø 15x 5mm dimensions were used. They were purchased from 
the firm ENES [22]. Their assembling consisted in uniform 
circumferential location of them in an appropriate number 14; 
10; 18 pieces, respectively. 


The basic features of the permanent magnets used in the 
tests are given below: 
- magnetic materialNeodymium N38 
- residual magnetic induction B, = 1.23 T 
- coercive force H, = 912 kA/m 
- maximum energy density (BH),,,, = 294 kJ/m3 
- maximum working temperature 150 °C 


Sealing of fermentation tanks, sealing at | For speakers of medium and 


high motion speeds, for food industry 


For speakers of low 


full range as well as sirens tone and sirens 


Pole shoes and sealing stage sleeves 


The pole shoes and sealing stage sleeves were manufactured 
of a low-carbon steel with the relative magnetic permeability 
u, > 2000. 


Magnetic fluid feeder and tool for the 
assembling/disassembling of sleeves with sealing 
stages 


For dosing the magnetic fluid the laboratory feeder 
LabMATE LM 100 of HTL firm, was used [23]. The feeder 
makes it possible to dose the fluid in the amount from 10 to 
100 ul with the filling accuracy error of + 0.8 % and repeatability 
error of < 0.3 %. 
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For the filling process the spouts of 200 ul capacity made 
by UNITIPS® firm were applied. 

To apply magnetic fluid to multi-stage seal it was necessary 
to use a special device as the sleeve with sealing stages was 
subjected, during its mounting, to action of a very strong 
magnetic field. The device was aimed at stabilization of the 
sleeve with sealing stages in a position which could make 
it possible to deliver magnetic fluid to successive sealing 
stages. 


PROCEDURES FOR CONDUCTING THE 
TESTS AT THE MFSL 135 TEST STAND 


To obtain repeatability of results of measurements aimed 
at determining critical pressure, critical motion speed and life 
of magnetic fluid seal it was necessary to elaborate and apply 
uniform testing procedures. Such procedures applicable to all 
the tests in question are described below. 


Procedure of test stand preparation to conducting 
the tests 


The procedure which makes it possible to conduct the tests 
consists in mounting the sleeve in an approporiate sequence, 
ends with application of a magnetic fluid and filling the testing 
chamber with utility water. 


Procedure of determination of critical pressure 
of the seal 


1. After 60 min passing from magnetic fluid application, 
to activate measurement mode at the control panel and 
Labview software, open the valve cutting-off compressed 
air inflow, set-up the flow throttling valve to the pressure 
increase rate within the range of (0.001+0.005) MPa/s. 

2. To increase the pressure till a leakage from the seal appears 
and simultaneously observe possible drop of pressure. 

3. After ending the test, to disactivate the measurement mode 
of LabView software. 

4. Toread, from the Labview records, the maximum pressure 

value observed during the test (i.e. the pressure at which 

penetration of the seal occurred). 

To introduce the measurement data into the testing card. 

To empty and open the testing chamber. 

To dismount and clean, by using extraction naphtha, all 

elements of the seal and testing head. 


Soo 


Procedure of determination of critical motion 
speed of the seal 


1. To activate measurement mode at the control panel and 
LabView software. 

2. To set up frequency converter to the maximum frequency 
value of 100 s-! and speed rise time of 250 s. 

3. After 45 min passing from magnetic fluid application, to set 
up water pressure in the chamber in compliance with the 
relative pressure value p/P = 0.8; 1.6; 2.4, respectively 
— in line with the testing programme. 

4. After 60 min passing from magnetic fluid application, to 
start up the water cooling circuit of the sealing system (for 
seal temperature stabilization), and next — start up the drive 
motor. 

5. To monitor increasing rotational speed of the motor and 
occurring leakage due to lack of tightness of the magnetic 
fluid seal. 
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6. After occurrence of lack of tightness to stop the motor and 
LabView software work. 

7. From Labview software records and monitored leakage 
occurrence instant, to determine the maximum value of 
rotational motion speed of the seal. The speed at which the 
leakage occurred is that critical for the seal. 

8. To introduce measurement data into the testing card. 

9. To empty and open the testing chamber. 

10.To dismount and clean, by using extraction naphtha, all 
elements of the seal and testing head. 


Procedure of determination of seal service life 


1. To activate measurement mode at the control panel and 
LabView software. 

2. To set up frequency converter to a given rotational speed 
in line with the testing programme, and speed rise time of 
20s. 

3. After 45 min passing from magnetic fluid application, to set 
up water pressure in the chamber in compliance with the 
relative pressure value p,/p,,, = 0.8; 1.6; 2.4, respectively 
— in line with the testing programme. 

4. After 60 min passing from magnetic fluid application, to 
start up the water cooling circuit of the sealing system (for 
seal temperature stabilization), and next — start up the drive 
motor. 

5. To monitor increasing rotational speed of the motor and 
occurring leakage due to lack of tightness of the magnetic 
fluid seal. 

6. After occurrence of lack of tightness to stop the motor and 
LabView software work. 

7. From Labview software records and monitored leakage 
occurrence instant, to determine the leakage occurrence 
time passing from the start-up of the test. The time passing 
from the start-up of the test to the first leakage instant is the 
seal service life in the conditions of the test in question. 

8. To introduce measurement data into the testing card. 

9. To empty and open the testing chamber. 

10. To dismount and clean, by using extraction naphtha, all 
elements of the seal and testing head. 


Measurement of pressure in the MFSL-135G 
testing head 


In the MFSL-135G testing head five pressure measuring 
channels are installed. Channel 1 is intended for the measuring 
of compressed air pressure. Channels 2 through 5 serve 
to measure water pressure in the testing chamber. To this 
end in the measurement disc are made two radial holes of 
different depth, connected to the testing chamber, that makes 
it possible to measure pressure values at different distances 
from shaft rotation axis. The measurements made by using 
Channel 2 are especially important for the performed tests as 
the outlet of its measurement hole is located directly in the 
neighbourhood of sealing stages filled with magnetic fluid. 
According to the programme realized in the frame of this 
project the pressure measurements from Channel | and 2 were 
used. Channel | serves to control pressure values set up in the 
phase of preparation of the test stand for carrying out the tests. 
Channel 2 makes it possible to measure water pressure just at 
the magnetic fluid seal. In Fig. 7 is shown an example run of 
the pressure in the testing head during the statical penetration 
test recorded by the LabView software. The record comes 
from Channels 2 through 5, in accordance with the pressure 
measurement points. 
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Fig. 7. An example run of the pressure in the testing head during the statical 
critical pressure test. The image taken from the graphical interface of 
LabView software, for four channels (Channels 2 through 5); the pressure 
measurement points — as above described 


SUMMARY 


In compliance with the method of conducting the tests it 
was assumed to determine, for each of the tested version of 
the seal, the following: 

- statical penetration pressure, 
- limiting working speed. 


Each of the measurements provides significant information 
which makes it possible to reach the assigned research aim, 
1.€.: 

- elaboration of a seal life curve and determination of the 
highest motion speed for which the seal is still capable of 
reaching, without any failure (lack of leakage), the basic 
number of cycles (rotations) equal to 1 mln, as assumed in 
these tests. 


A very important factor for the performed tests is the water 
pressure set-up in the measuring head. The factor was taken into 
account in the form of the relative pressure defined as the ratio 
of the set-up pressure and the critical pressure, (p,, = PPn, 
defined as the seal penetration pressure per one sealing stage, 
obtained in the statical test. 

Due to the assumption of such concept of setting-up testing 
pressure of water it was possible to simply take into account 
influence of the factor on seal life for particular tested magnetic 
fluids, however the assumption also makes comparison 
of service life of different fluids difficult because of large 
differences in values of the set-up statical testing pressure for the 
tested magnetic fluids, which result from very large differences 
in values of their statical penetration pressure. For instance, the 
highest long-term operation speed equal to 45 rps was obtained 
during the tests with the FLA 002.25 magnetic fluid at the set- 
up pressure of 0.0237 MPa (p,, = P/P = 9.8), resulting from 
the low penetration pressure of the fluid equal to 0.0297 MPa 
per one sealing stage. Since the analogous test at the value 
Py = P/Pxu = 9.8 for the FLA 003.045 fluid requires to set-up 
the water pressure value of 0.058 MPa (the penetration pressure 
per one sealing stage for the fluid is equal to 0.072 MPa) then 


the long-term speed for the seal filled with the fluid amounts 
only to 20 rps. The results indicate a significant influence of 
pressure on service life of magnetic fluid seal operating in water. 
Limiting working speed of the seal brings significant research 
information as it makes it possible to assess the speed value 
from which elaboration of seal life curve should be commenced, 
and consequently, determination of seal life for the basic 
number of cycles equal to 1 mln. However from the utility 
point of view the parameter is rather not important. 

The tests conducted with different numbers of magnets in 
the sealing system has to show influence of size of magnetic 
field on limiting speed of magnetic fluids. 

The preliminary results of seal life tests for majority of 
the tested magnetic fluids and seals provide similar picture 
on the diagrams. At higher rotational speeds short working 
time to leakage (low number of cycles) is observed. At lower 
motion speeds the time to leakage appears much longer and 
clearly tends to a speed value at which seal operation can be 
considered long-term. The full results of series of the tests will 
be described in the other part of the article. 
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The effect of liquid temperature and viscosity 
on the amplitude-frequency characteristics of 
a viscotic torsion damper 


Wojciech Homik, Ph. D. 
Rzeszów University of Technology, Poland 


ABSTRACT 


The article discusses causes of the appearance of transverse, longitudinal and torsional 
crankshaft vibrations in multi-cylinder internal combustion engines. Particular attention 
is paid to the torsional vibration which is the most severe threat to engine crankshafts. 
Damping methods making use of torsion dampers are presented. With the reference to 
viscotic dampers, problems with their damping efficiency are discussed in the context of 
viscosity changes of the damping liquid. The article also presents the amplitude-frequency 
characteristics of a series of viscotic dampers, which were recorded experimentally on the 
research rig and on a real object. An idea of vibration damper metric is given. 


Key words: torsional vibration; torsional vibration damping; torsional vibration dampers; 
crankshaft; amplitude-frequency characteristics; liquid viscosity; damper metric 


INTRODUCTION 


The internal combustion engine. Could the designers of the 
first internal combustion engines predict that their construction 
would become the source of progress and development of 
transport facilities? 


a "Z Ra 
> 4 
X <M 
Fig. 1. Designers of piston internal combustion engines [12]: a) Nicolaus- 
August-Otto, b) Eugenem Langenem, c) Etienne Lenoir, d) Rudolf Diesel 


Surely, none of them expected that the basic idea of their Fig. 2. Design of Otto engine of 1876 
construction (Fig. 2, 3, 4) [12] which changes the reciprocating - ii se 
linear motion into the rotary motion would remain in ; 
a practically unchanged form until now. 

The forces acting on the engine crankshaft system 
include: 

e gas pressure forces P,, 
e inertial forces P,. 


The gas pressure forces P, and the inertial forces P, are 
periodically changing forces. In four-stroke engines the period es J 
of the inertial forces is equal to one full revolution of the ii A FS 
engine crankshaft, while that of the pressure forces — to two Fig. 3. Construction of Etienne Lenoir engine 
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Fig. 4. Construction of Rudolf Diesel engine 


revolutions. In two-stroke engines the latter period is also equal 
to one full crankshaft revolution. 

The resultant force of the above forces is their sum [1, 13, 
14, 15, 16] 


P=P,+ P, (1) 


The analysis of changes of the resultant force P as a function 
of the crank angle @(t) reveals that at the final compression 
phase and the initial expansion phase the inertial forces unload 
the crankshaft-piston system, while in the remaining time they 
are the sources of extremely high load, especially during the 
suction and exhaust strokes. 

Since in general, the axis of the connecting rod is not in line 
with the cylinder axis during engine operation, the resultant 
force P acting on the piston is not transmitted solely by the 
connecting rod. It is shifted towards the piston pin axis, in 
which it is decomposed into two components: 

e N component acting on the cylinder surface, 
e S component acting along the connecting rod axis. 


The magnitude of the force N depends on the coefficient À, 
being the ratio of the crank arm length r and the length L of the 
connecting rod (the larger A, the larger the force N). 

The force S acting along the connecting rod axis acts on the 
crank, where it is again decomposed into two component forces: 
the T component which is tangential to the circle described 
by the crank, and the radial component K acting along the 
temporary crank arm position. 

The above forces are given by the relations: 


P 
on sin(ọ+ß) (2) 


T=Ssin(9+)= 


K =Scos(9+f) = 


= soD © 


cos 
where: 
@ — crank angle, 
B — angle of crankshaft deflection from cylinder axis. 


72 POLISH MARITIME RESEARCH, No 4/2012 


The periodic nature of changes of the gas pressure forces 
P, and the inertial forces P, is the source of the following types 
of crankshaft vibrations [1, 14, 15]: 
e transverse vibration, 
e longitudinal vibration, 
e torsional vibration. 


These vibrations are sort of protection measure, taken by the 
engine components made of elastic materials, which consists 
in yielding to the applied load by gradual absorption of the 
transmitted energy in the form of vibrations. Many engine 
components, including the abovementioned crankshafts, would 
be damaged very fast if they did not have the ability to absorb 
the energy in the form of elastic deformations (Fig. 5). 


Crank at rest 


Deformation caused 
by bending 


Deformation caused 
by torsion 


Fig. 5. Sample elastic crank deformations caused by T force 


VIBRATIONS 


Transverse vibration 


The transverse vibration of the crankshaft, also sometimes 
referred to as the bending vibration, is caused by the action of 
periodically changing radial forces K on the crank pins. This 
vibration results in bending the shaft between the bearings 
(Fig. 6). 


Fig. 6. Shaft deflection between bearings 


Longitudinal vibration 


The longitudinal vibration of the crankshaft is directly 
related to its transverse vibration. Any shaft deflection causes 
its axial shift (Fig. 7). 

In most cases this vibration does not disturb the engine 
operation (especially in motor-car multi-cylinder internal- 
combustion engines) and is not a severe threat to the durability 
of the crankshaft, which reveals relatively high rigidity. 


Fig. 7. Longitudinal crankshaft vibration 


Piston 


Crankshaft Thrust bearing 


Screw 
propeller 


Line shafting 


Fig. 8. Real and virtual propulsion system of a watercraft 


It should be stressed, however, that this vibration can be 
a severe problem in marine large-power engines. They make 
the entire structure (Fig. 8), consisting of the engine crankshaft, 
the flywheel, the line shafting and the screw propeller, move 
periodically along its axis. The amplitude of the longitudinal 
vibration of this structure depends in practice on individual 
design solutions of the screw propeller, in particular on the 
number of paddles [14]. 


Torsional vibration 


Among all forces acting on the piston-crankshaft system, 
the force which makes the crankshaft rotate is the T force 
tangential to the circle described by the crank. 

This force affects the engine torque which is equal to 


MT r= rsin A 
cosB 


where: 

M(@) — shaft torque, 
T(@) — tangential force, 
r — crank arm length. 


The T force, which is one of two components of the S force 
acting along the connecting rod axis, is a periodically changing 
force. Its period is equal to 2r for two-stroke engines and 4r 
for four-stroke engines. The changing T force is a source of 
accelerations in the rotational motion of the engine crankshaft, 
and torsional vibrations, as a further consequence (Fig. 9), which 
change with the changed rotational speed of the shaft [15]. 

Among all vibrations, the torsional vibration creates the 
biggest threat to the crankshaft. 

It is noteworthy that the crankshaft deflection is only limited 
by the torsional stiffness of the shaft, and the amplitude of the 


torsional vibration can exceed the permissible values. In case 
of total lack of damping, the curve according to which the 
vibration changes tends to infinity for any speed generating 
a resonance harmonic. 


Fig. 9. Torsional crankshaft vibration 


Torsion break of the shaft takes place when the permissible 
value of the amplitude @,,, (permissible torsion angle) is 
exceeded 


Piz > Ddop (5 ) 

The permissible value of the torsion angle @,,, depends 

on the conditions of engine operation and is determined by 

general tolerance of geometrical parameters. A term which 

is frequently used for its easy comparison is the elementary 
torsion angle [2, 3]: 


1 GIhlm 
where: 
M, — torque, 
G — torsional modulus, 
I, — polar moment of inertia 
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For steel propeller shafts and fluctuating loads: 


top = 0 004) =£ | (6a) 
m 
while for alternating loads 
P'iop = 0.0025 H (6b) 
m 


The permissible values of the torsion angle are most 
frequently within the limits: 


P'iop = (0.002 +0. on] =) (6c) 


For the shaft revealing changing torsional stiffness (profiled 
crankshaft) the basic condition (6) takes the form: 


oe Oop (7) 


= 


TYPES OF TORSION DAMPERS 


Detecting the torsional vibration of the engine crankshaft is 
more difficult than in case of other vibration types. Its presence 
can be sometimes testified by, for instance, non-uniform engine 
operation, the timing gear of which is driven by the torsionally 
vibrating crankshaft via a mechanical transmission (for instance 
belt, chain, or gear transmission). 

The variety of vibration forms and the existence of different 
harmonics of the tangential force T are the reasons why the 
crankshaft can work in the zone of resonance at different 
rotational speeds of the engine. 

The engine operation in critical (resonance) conditions 
is most frequently avoided by the use of vibration dampers 
(Fig. 10). 


Fig. 10. Torsion dampers 


Their task is to decrease the amplitude of the resonant 
torsional vibrations of the engine crankshaft. A properly 
designed (selected, “tuned”) torsion damper provides 
opportunities for decreasing the amplitude of torsion vibration 
resonance by even as much as ten times. But we should keep 
in mind that each damper consumes part of the effective power 
of the engine. 

In the past, the threat coming from torsional vibrations of 
watercraft power transmission systems was minimised using 
the following types of dampers: 

e frictional dampers, 
e rubber dampers, 

e viscotic dampers, 
e spring dampers. 
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The above named dampers are typical dynamic resonant 
vibration dampers in which the vibration is damped using 
the inertia force [4, 5]. Despite a common name, the dynamic 
dampers differ between each other not only by constructional 
designs but, what is more important, also by operating 
characteristics [10]. It is also noteworthy that these dampers 
work effectively only within a certain range of rotational 
speed. 

Among all the above mentioned dampers, it is the viscotic 
dampers which have found the widest application in watercraft 
power transmission systems. This most likely results from 
their simple structure, operating reliability, relatively simple 
manufacturing technology and low production costs, as 
compared to other solutions. 

Nowadays the viscotic torsion dampers are filled with 
stabilised silicone oils of AK type produced by Bayer and 
Wacker (Fig. 11). A basic parameter characterising the oil (like 
any other real liquid) is its viscosity, defined as the internal 
resistance which appears during displacement of a liquid layer. 
In this sense the viscosity is decisive about fluid friction: the 
higher viscosity, the larger the internal friction of the liquid. 


Fig. 11. Process of filling the viscotic torsion damper [8] 


The silicone AK oils are watery clear liquids, the viscosity 
of which changes with shear velocity changes (Fig. 12). 
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Fig. 12. Kinematic viscosity logarithm vs. shear velocity logarithm [20] 


Changes of liquid viscosity are caused not only by the 
increase of its shear velocity but also by changes of the 
operating temperature and oil ageing processes [6, 7]. Cases are 
well-known in which improper use of the damper has resulted 


in remarkable increase of the oil viscosity and change of the 
oil aggregation state (Fig. 13) [9, 10] 
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Fig. 13. Samples of oils — magnified 200 times: a) oil with particles of 

damper parts, b) oil having the consistence of a solid body 


Oil viscosity changes cannot be neglected at the stage of 
selection (designing) of the viscotic damper for a certain power 
transmission unit [11]. When selecting the damper, calculations 
are performed to determine an optimal liquid viscosity, i.e. the 
viscosity for which the value of the energy dissipated by the 
damper in given operating conditions reaches the maximum. 

As a rule, the damper calculations take into account the 
nonlinearity of the liquid. Among other publications, this issue 
was discussed in detail in the work entitled “Optimisation of 
the viscotic torsion damper” [4]. In that work, in order to define 
the effect of liquid nonlinearity of the efficiency of the viscotic 
torsion damper, the damping coefficient a which, as mentioned 
above, depends on the liquid viscosity, was described using the 
following functions: 


= An l-ci. -Pp Oy €< 0,00) (8) 


O=Ay expl- clp, -Pp 


), c20 (9) 


=O yo (10) 
where: 
Oy, — damping coefficient of a damper with Newtonian 
liquid, 
Po — Mp — relative speed of the casing with respect to the ring, 
c — coefficient taking into account the effect of the shear 


velocity on viscosity. 


Sample results of calculations are shown in Figs. 14 - 18. 
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Fig. 14. Casing vibration amplitude vs. logarithm c 
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At present, for technological reasons the viscotic torsion 
dampers make use of oils the viscosity of which does not 
exceed 300000 [cSt] at temperature of 20°C. In those cases, 
when modelling the viscotic torsion dampers we can assume 
that the analysed liquid is the linear liquid. 
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Fig. 15. Casing vibration amplitude as a function of ay. 1: € = 0; 2: c = 
10-6; 3: c 
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Fig. 16. Amplitude-frequency characteristics for different c values 
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Fig. 17. Work L vs. kinematic viscosity 


Before the damper is finally allowed to be used in 
a given power transmission unit, each respected producer of 
dampers looks for confirmation of his theoretical analyses in 
experimental tests, also those performed on a real object, if 
possible. The results of these tests provide opportunities for 
final verification of the already performed theoretical analyses 
of vibration dampers. 

The tests are oriented on determining the amplitude- 
frequency characteristics of the examined dampers. These 
characteristics, complemented by such damper parameters as: 
e I ,— mass moment of inertia of the ring in the damper, 

e I,,— mass moment of inertia of the casing, 
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e 1 - viscosity of the silicone oil, 
compose the “metric” of the damper, which should be made 
available by each producer of vibration dampers. 
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Fig. 18. Damper efficiency as a function of kinematic viscosity v for liquid 
with n = 0.9 


Damper tests which make it possible to work out their 
amplitude-frequency characteristics are most frequently done 
in temperatures close to those in which the damper is expected 
to work most often (different operating conditions are expected 
for a damper used in watercraft power transmission systems or 
electric current generation units, and different for that working 
in power transmission systems in vehicles). 

In order to illustrate the differences in amplitude-frequency 
characteristics between different dampers, a series of viscotic 
dampers: 

e TW 30 -viscotic damper with damping liquid of viscosity 
of 30000 [cSt], 

e TW 60- viscotic damper with damping liquid of viscosity 
of 60000 [cSt], 

e TW 100-viscotic damper with damping liquid of viscosity 
of 100000 [cSt], 

e TW 130-—viscotic damper with damping liquid of viscosity 
of 130000 [cSt], 

e TW 160-viscotic damper with damping liquid of viscosity 
of 160000 [cSt], 

e TW 200-viscotic damper with damping liquid of viscosity 
of 200000 [cSt], 

e TW seized, 

having the same geometric parameters were tested at three 

different operating temperatures: 

e ambient temperature 25°C, 

* mean operating temperature 45°C, 

e saturation temperature 80°C, 

within the range of the exciting frequencies changing between 

25 and 300 [Hz]. These frequencies were introduced to the shaft 

of the research rig using the motor TYP FD. The parameter 

which was measured in the tests was the torsional vibration of 

the damper hub which was directly connected to the vibrating 

shaft. 

A brief analysis of the presented results reveals that for each 
damper at a given temperature there is an optimal viscosity 
value for which the shaft torsion angle reaches it minimum. It is 
noteworthy that within the range of lower exciting frequencies, 
between 25 and 100 [Hz], the viscosity of the damping liquid 
does not affect remarkably the efficiency of the damper, 
irrespective of the temperature. 

The viscosity of the liquid is decisive for the efficiency of 
the damper within higher excitation frequency ranges, between 
200 and 300 [Hz]. 
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Fig. 19. Comparing amplitude-frequency characteristics of the TW200 
damper in temperatures 25°, 45°, 80° [C] 


The amplitude-frequency characteristics provide good 
opportunities for final selection of optimal viscosity of the 
liquid in a damper which damps torsional crankshaft vibrations 
in an electric current generator unit which most often works 
with constant rotational speed (and practically constant 
excitation frequency). It is much more difficult to select the 
optimal liquid viscosity for a damper used for damping torsional 
crankshaft vibrations in a marine engine and, perhaps most 
of all, a motor-car engine in which the excitation frequency 
changes with rotational speed changes. In those cases the 
practically optimal viscosity is calculated for most dangerous 
excitation frequency (resonance). 

The optimal damping is calculated using the formula: 


Qop = 1,0, (7.4) 


where: 

real frequency at which hub vibration is the smallest, 
mass moment of inertia of the inertial ring. 
For instance, the optimal viscotic damping @, for the 
damper TW130 at temperature of 45°C is: 


(45) = 0.006750877:1650 = 11.195 [Nms] 


opt 


Ooptrw130 
CONCLUSIONS 


- In order to confirm the theoretical analyses and the result 
obtained in the experiments done on a model research rig, 
investigations were performed on a real object, which was 
a car VW Passat with the engine cubic capacity of 3.6 | 
and power of 290 [KM]. A parameter which was recorded 
during the tests performed in the engine test house was the 
amplitude of the torsional vibration of the free shaft end, 
within the engine rotational speed ranging between 700 and 
6100 [rev/min]. 

- The harmonic analysis of the signal recorded during the 
tests has revealed, most of all, that the examined dynamic 
system cannot safely work without a torsional vibration 
damper, especially when the rotational speed exceeds 3500 
[rev/min], and that the biggest threats for the system are the 
3-rd and 6-th harmonics. 

- Moreover, the tests performed on the real object have 
confirmed the conclusions resulting from the theoretical 
analyses and research rig tests about the existence of 
optimal values of liquid viscosity [4] for which the maximal 
resonance amplitudes of the engine crankshaft reach their 
minimum (Fig. 20, 21). 
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Fig. 20. Shaft torsion angle changes vs. oil viscosity 
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Fig. 21. Shaft vibration frequency changes vs. oil viscosity 


- The analysis of the characteristics shown in Figs. 20 and 
21 clearly confirms the opinion that the classical viscotic 
vibration dampers are the optimal dampers only in a very 
narrow range of frequency excitations, which is caused by 
the fact that during the operation of the dynamic system we 
cannot change physical parameters in a viscotic damper (its 
physical characteristics are constant in practice and cannot 
be adapted to changing working conditions — changing 
excitations). That narrow band of efficient operation of 
these dampers is undoubtedly their disadvantage, but their 
simple structure and relatively low production costs are the 
reasons why they are commonly used for damping torsional 
crankshaft vibrations in power transmission systems. 
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Data fusion in a navigational decision support 
system on a sea-going vessel 
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ABSTRACT 


The problem of data fusion in a navigational decision support system on a sea-going vessel 

has been analyzed. The computing algorithm herein applied for solving the formulated 

problem is based on a multi-sensor Kalman filter. On the practical side, results of the 

tests done in real conditions are shown. The tests conducted onboard m/s Nawigator XXI, 

have been aimed at the verification of the proposed computing algorithm implemented in 
a prototype navigational decision support system. 


Key words: data fusion; navigational decision support system; multi-sensor Kalman filter 


INTRODUCTION 


In the literature on the problems of information processing 
the terms data integration and data fusion are very often used 
as synonyms [9]. This author assumes that the term data 
integration is understood as a process in which data from 
various, generally physically separate sources are combined 
into a consistent and readable whole in a manner allowing to 
obtain new, more comprehensive information that would be 
impossible to get from each individual source. Data fusion, 
on the other hand, is understood as a process of merging data 
from various sources, however, describing the same quantity. 
Therefore, this process results in obtaining more reliable and 
more accurate information. 

Data fusion algorithms are naturally applied when input 
data needed to solve a specific problem are simultaneously 
obtained from many redundant measurement devices. For 
further overall processing these data have to be transformed into 
one signal, which will take account for the information from all 
sensors. Thus obtained data should be more precise than those 
derived individually from one source. Besides, this approach 
increases the level of safety of system which may be affected 
by a possible lack of input data. Even if one source (device) 
fails the others guarantee effective operation. This observation 
is of particular importance in real time systems. 

Methods and algorithms of data fusion find applications 
in many problems across various scientific disciplines, e.g. 
in medicine [7], land transport [14], processes of knowledge 
retrieval from data bases [17], problems of mobile agents [18], 
satellite image integration [11]. In this article data fusion is used 
for solving problems of dynamic object (vessel) positioning in 
marine navigation. 
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Chapter V, paragraph 19 of the SOLAS (Safety Of Life 

At Sea) Convention [5] deals with mandatory fitting of sea- 
going ships with navigational systems and equipment. Which 
equipment items a ship should carry depends on the date the 
ship’s building commenced (laying the keel) and on its gross 
tonnage. According to the above mentioned document, ships 
with a tonnage above 500 GT (majority of sea-going vessels) 
whose construction started after 01 July 2002 should be 
equipped with, but not limited to: 

- Global Navigation Satellite System (GNSS) receiver, 
determining ship’s position (latitude and longitude) 
automatically, 

- AIS system (Universal Shipborne Automatic Identification 
System), providing, inter alia, information on own ship’s 
position (latitude and longitude). 


Practice has shown that most sea-going vessels comply with 
the above requirements as they have redundant equipment. This 
is an effect of continually reduced costs of GNSS receivers 
due to growing demand for satellite navigation systems in 
land transport. 

No wonder, modern ships have several different devices 
with sensors measuring their own position (latitude and 
longitude). The navigator, having differing data on his ship’s 
position, may face difficulties in making the right decision. 
On the other hand, relying exclusively on measurements 
from one device is a risky solution as the information source 
may fail. Therefore, it seems purposeful to apply the process 
of navigational data fusion. One of the methods to solve 
this problem is application of a multi-sensor Kalman filter 
[1, 15,16]. 


ALGORITHM OF NAVIGATIONAL DATA 
FUSION WITHIN THE KALMAN FILTER 


To meet the demands of the herein considered algorithm 
of navigational data fusion, we introduce the following 
mathematical model of discrete stochastic system with | sensors 
[2,3,4]: 

x(tl)+ = @-x(t)+w(t) 


y(t) =H, -x(t)+v,(t) 1=1,2,K 1 a 
where: 
x(t) 
x(t) = a eR’ — state vector at instant t, 
v(t) 
(x(t), y(t)) — Cartesian coordinates of vessel position 


at instant t, 


(v,(t), v0) 


y(t) € R? 


1010 
[0101 
10010 

0001 


wit), vi(t) 


— vector of vessel velocity in the reference 
system at instant t, 

—measurement vector of i-th sensor 
(Cartesian coordinates of vessel position 
from a GNSS receiver, vessel velocity 
vector from a Doppler log) at instant t, 


— constant system matrix, 


— constant matrix of i-th sensor (GNSS 


. 1000 
receiver: i 0 | Doppler log: boo °p. 


— vectors of disturbances with characteristics 
of white Gaussian noise with zero 
expected values and covariance matrices 
Q and R,, respectively, at instant t. 


GNSS, GNSS, ... GNSS, LOG, 


GNSS awsstl 
B, (t) Bs rss (t) Viona (t) Vion (t ) 
L, (t) Lia (t) Viet mest (t) Viat (t) 


transformation transformation transformation 
geographical c. || geographical c. sss geographical c. 
Cartesian c. Cartesian c. Cartesian c. 


transfering to 
` | reference system 


y,(t)= Finn (t)= Visit (t)= 


BEAD ib) Vran t) a (r) 
= H j j = o 7 Pa (t) i Va (t) 


Kalman filter Kalman filter | 
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weighted 
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navigational decision 
support system 


Fig. 1. A general diagram of the proposed navigational data fusion algorithm 
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A general diagram of navigational data fusion algorithm within the Kalman filter is shown in Fig. 1. 
Denotations of the quantities included in Fig. 1: 
lonss — number of GNSS receivers (1 < lensg < 1), 
B,(t), L(t) — latitude and longitude of the vessel obtained by i-th GNSS receiver at instant t (1 < i < lenss) [°]; 
Vion,(t), Vion(t) — longitudinal and transverse velocities of the vessel obtained by i-th Doppler log at instant t (Igyss + 1 < i <1) 
[m/s], 
(x; (t), y,(t)) — results of measurement of the Cartesian coordinates of the vessel, obtained by i-th GNSS receiver at instant t 
and transferred to the position of the first GNSS receiver (2 <1 < Igngg) [m], 
V(t), Vy(t)  — results of measurement of vessel velocity vector obtained by i-th Doppler log at instant t (Igyss + 1 < i < 1) 
[m/s], 
(0) — estimates of state vector obtained through Kalman filtration for i-th subsystem, 
x(t) — vector of state estimates fusion. 
The following navigational systems and equipment are sources of input data for the algorithm, where they undergo the 
fusion process: 
- satellite navigation system receivers and the AIS system which carry out measurements of own vessel’s position (latitude 
and longitude), represented as GNSS receivers in Fig. 1, 
- Doppler logs, which provide measurements of longitudinal and transverse velocity components of own vessel. 


For the effective operating of the algorithm the vessel must carry at least two of the mentioned devices. Additionally, 
a gyrocompass allowing to measure own vessel’s heading needed in the process of internal data transformations, is required. 
The circulation of signals in the proposed algorithm (Fig. 1) will be further described in detail. 


Transformation of geographical coordinates into Cartesian coordinates 


Results of measurements of own vessel position delivered by GNSS receivers have a form of geographical coordinates 
(Fig. 2), so that they cannot be directly used in the model (1). To transform them into the Cartesian coordinates the Gauss-Kriiger 
projection was applied [6]. Then the relevant transformation formulas obtain the following form: 


x(t)= N(t)-cos(B(t))- AL(t)+ 
+ -N(t)-cos*B(t)-d-T’ E+ O ALOY + 


a N(t)-cos°(B(t))- 6-18 -T +T +14- 4) -58 T O-n 0) (ALAY + 


30 NG): cos B(t)-61-479-T? ()+179-T!)- T°) QL)’ 


y(t)= SO+5 — - NG): co? BOY- (AL@Y + 
E i 


TO 


“720 
TO 


40320 


(2) 


-N(t): cost BEY- 6= T?(t)+ 9-17 (+4 n O ALO + 
-N(t): cos°(B(t)-(61—58- T? + T4(t) + 270-1 (t)-330-T? (t)- n O): (ALE + 


+ —— 


-N(t)-cos®(B(t)) 0385 —3111- T?(t) +543- T ©- Tf €X- (ALOY 
where: 
N(t)= 


2 
a 
re 7 


y= ~ -cos*(B(t) = 


radius of the first vertical curvature (for WGS 84 ellipsoid we can assume: a = 6378137 [m], 
b = 6356752.3141 [m]), 


non-dimensional auxiliary variable, 


AL(t) = i -L, — difference between measurement results of a given longitude and the longitude of the axial 
meridian of the examined area (for the Szczecin area we can assume L, = 0.2618 [rad]), 

T(t) = tg(B(t)) — non-dimensional auxiliary variable, 

S(t) — length of meridian arc B(t) [m]. 
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Fig. 2. Data from a GNSS receiver 
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The meridian arc length B(t) is found from the following 
formula: 


S(t) = 6367449 - B(t) + 
— 16038.50874 - sin (2 - B(t)) + 
+ sin (4 - B) — 0.022198 - sin (6 - B(t)) + 
+ 0.00003 - sin (8 - B(t)) 


The proposed transformation reproduces accurately only 
measurements close to the axial meridian. The difference in 
the transformed length grows when it is further from the axial 
meridian, reaching a maximum at the edge of the zone (for 
a six-degree zone it equals 67 cm per | km of length). The area 
subject to navigational situation analysis to which the data from 
the fusion process will apply, should not exceed a circle with the 
radius of 10 nautical miles. If we take this into account and the 
fact that the selected axial meridian has to cross the analyzed 
area we can assume that transformation errors are minor and 
will not affect the quality of the fusion. 


(3) 


Transferring the measurement results from 
GNSS receivers to one common point 


Aerials of GNSS receivers installed on board a ship are 
generally mounted in different places. As the aerials are some 
distance apart, the obtained measurement results should be 
brought, after passing through the block of transformation 
from geographical coordinates to Cartesian coordinates, to one 
common place to be subjected to the fusion process. This is due 
to the fact that in the model (1) the vessel position is understood 
as a specific material point. Most often it is interpreted as the 
hull centre of gravity. Unfortunately, it rarely happens that any 
of the aerials is placed exactly at such point, and the process of 
bringing the measurement results to a specific point is always 
burdened with an error (this is due to the need of using, in the 
process, the gyrocompass heading which is a measured quantity 
burdened this way with an error). Therefore, it is decided that 
in the proposed algorithm the place to which all measurement 
results are transferred is the position of one GNSS receiver 
aerial, preferably that of the most accurate receiver (if there is 
such). In this approach measurement results from one GNSS 
receiver (GNSS,) will be free from additional deformations 
caused by the geometrical transformation. This transformation 
consists in combining the rotation by the angle (course) 
and the translation by the vector according to the following 
equations: 


X(t) = x(t)+ w, cos(y(t))— w, sin(y(t)) 
¥(t)= y(t)+ wy cos(y(t))+ wy sin(y(t)) 


where: 

w = [Wp w,] — vector which specifies the position ofa GNSS 
aerial relative to the GNSS, receiver aerial 
(Fig. 3) [m], 

y — gyrocompass heading [rad]. 


Shifting the measurement results obtained from GNSS 
receivers will increase the measurement error by the value: 


2 2 ow Cg ro 5 
2A W +W, snf A (5) 
where: 


e — gyrocompass error [rad]. 


gyro 


As the gyrocompass error should not exceed 0.009 rad, 
while the distance between the GNSS aerials should not 


exceed 15 m, we can assume that errors caused by bringing 
the measurement results to one common location are negligibly 
small and will not affect the fusion quality. 


Transferring the velocity measurement results to 
the reference system 


In parallel with position measurement transformations, 
similar operations concerning velocity (over ground) 
measurements should be performed by using Doppler logs. 
Measurements of longitudinal and transverse velocity 
components of a vessel should be related to the introduced 
reference system through a rotation by the gyrocompass 
heading in accordance with to the following equations: 


Vy (t) = Vion cos(y (t »)) —Viat sin (y W) 
Vy (t)= Viat COSQY(E)) + Vion Sin Qy (t)) 


The transformation will increase the measurement error 
by the value: 


_ fe 
U + ve sin( “| (7) 


Since the gyrocompass error should not exceed 0.009 rad 
and the vessel speed should not be higher than 15 m/s, we can 
assume that errors caused by bringing velocity measurements 
to the reference system are also negligibly small and will not 
affect the fusion quality. 


The Kalman filter 


Estimates of the state vector %,(t) for i-th subsystem 
(relating to a sensor) are obtained by using a Kalman filter [8] 
for the model (1): 


&,(t\t-1)= ®-&, t-1\t-1) 

EO =y -H-4 t\t-1) 

P.(t\t-1)= ®-P.(t-1\t-1)-®'+Q (8) 
K.(t) =P (t\t-1) HD [H P (t-11) HR] 
P(t\t) = P(t\t-1)-K,(t)-H, -P,(t\t-1) 

XK (t\t) =x(t\t-l+K,(t)-E,® 


(6) 


where: 
E(t) — vector of innovation at instant t, 
K((t) — filter gain matrix at instant t, 


P(t\t— 1) — filtration errors covariance matrix determined 
without the knowledge of filter gain at instant t, 

P (tt) — updated filtration errors covariance matrix, 

x,(t\t — 1) — estimate of the state vector determined without the 
knowledge of measurements at instant t, 

X,(t\t) — updated estimate of the state vector. 
Data fusion 


The fusion of a set of data obtained from 1 sensors is 
expressed as a weighted mean: 


X(t) =A,(t)-&,(t)+ A,(t)-X,(O+K+A,(1)-% (1) O 


where: 
A(t) — weight matrices. 


The weight matrices are derived from the formula: 


A= LPO PO (10) 
jal 
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where: 
P(t) | — matrix of cross covariance of error filtration between 


i-th and j-th subsystem of model (1). 


Matrices of the cross covariance of filtration errors are 
determined from the following formula: 


P,(t)=[1, -K,(0- Hj] -[® -P, (t-1)-©'+ Q]} 


T (11) 
1, -K;6-B;] 
where: 
L — 4x4 unit matrix. 

The proposed approach is optimal as it minimizes the trace 
of error variance matrix of the fusion estimator [15, 16]. 

It should be emphasized that due to the used model (1), the 
algorithm of navigational data fusion is a positioning algorithm 
for rectilinear movement. That is why the fusion process cannot 
be continued when the vessel starts turning. 


DATA FUSION IN THE NAVIGATIONAL 
DECISION SUPPORT SYSTEM 


Rapid advancements in information technologies we 
have been witnessing in recent years, significantly affect the 
organization and execution of maritime transport processes. 
More attention is being paid to the construction of intelligent 
systems which utilize a wide range of telecommunication, 
computer, automation and measurement technologies, in order 
to protect traffic participants, increase effectiveness of the 
maritime transport system and protect natural environment 
resources. Intelligent transport systems (ITS) [12] are advanced 
applications with a capability of safer and better coordinated 
use of transport networks. One example of such applications 
is the Vessel Traffic Monitoring and Information System 
(VTMIS) presently being developed within the EU. The 
VTMIS consists of several components: the Vessel Traffic 
Service (VTS), Automatic Identification System (AIS), Ship 
Reporting System (SRS), Maritime Assistance Services 
(MAS), Long Range Identification and Tracking System 
(LRITS), computer-aided system of information exchange 
(SafeSeaNet). 

In December 2005 the International Maritime Organization 
(IMO) started working on developing and implementing a new 
navigational strategy - e-navigation [10]. E-navigation is 
understood as gathering, integration, transmission, management 
and presentation of navigational information by means of 
electronic formats for the navigational support of port-to-port 
operations. The concept of e-navigation goes in line with the 
concept of maritime intelligent transport systems. 

The development of navigational equipment and systems 
installed on sea-going ships created a demand for data integration 
and devising methods of navigational data presentation which 
would facilitate the navigator in decision making. Another stage 
consists in building specialized information systems assisting 
decision processes on the ship, such as pilotage and docking 
systems. Their development proceeds towards decision support 
systems which in turn have such functions as current situation 
analysis and generation of proposed solutions to dangerous 
(i.e.collision) situations. 

The research team at the Institute of Marine Navigation, 
Maritime University of Szczecin, has developed a prototype 
navigation decision support in open sea areas, thus making 
a contribution toward the global work on intelligent systems 
and meeting today’s demands of the maritime industry [13]. The 
computer-based system (hardware and software) can be easily 
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installed on board a ship and work in real time. The prototype 
has the following functions: 

- recording and decoding navigational data, 

- navigational data integration, 

- navigational data fusion, 

- identification of a navigational situation, 

- analysis and assessment of a navigational situation, 

- generating anti-collision manoeuvre, 

- prediction of a navigational situation, 

- presentation of the processed information. 


One of the functions is the process of navigational data 
fusion. In the prototype system this process is based on the 
implemented algorithm such as that presented in the previous 
chapter. The use of the computing algorithm ensures a number 
of advantages: 

- increased reliability of the system, 

- reduction of measurement errors, 

- continuity of the system operation, 

- data acquisition at greater frequency. 


THE RESULTS OF FIELD TESTS 


The algorithm of navigational data fusion implemented in 
the navigational decision support system was verified in real 
conditions on board the m/s Nawigator XXI. The field tests were 
carried out in the area of the Szezecin-Swinoujscie fairway 
at the southwest Baltic Sea. During the tests the ship had the 
following navigational equipment and systems essential for 
the algorithm verification: 

- AIS system - Nauticast X-Pack DS, 

- GNSS receiver - CSI MiniMax, 

- GNSS receiver - Koden KGP-913D, 

- GNSS receiver - Trimble NT 200D, 

- gyrocompass - Gyro STD22 Anschutz. 


For technical reasons the input data for the tested algorithm 
were signals from two GNSS receivers (MiniMax, Koden) and 
the gyrocompass. The MiniMax was chosen as the primary 
GNSS receiver (GNSS,) due to its higher accuracy. 

The disturbance covariance matrix was assumed as: 


0 0 0 0 

Q= 0 0 0 0 d 
0 0 0 0 
0 0 0 0 


while the sensor-related data (accounting for the characteristics 
of measuring devices) were assumed as: 


0.25 0 
i=l @ gre y 

225 0 a 
„= g a w, =[0.56 -1.04| 


These are the initial matrices of filtration error covariance 
and matrices of cross covariance of filtration errors (applied 
after each filter cancellation — subsequent start-up of the 
algorithm): 


10 0 0 
0 1 0 0 
P= (14) 
0 0 0.0625 0 
0 0 0 0.0625 


First measurements from each device were assumed as the 
remaining initial values. 

Fig. 4 and 5 present example records of ship positions 
measured when the ship proceeded along a straight track 
section. It can be seen that the measurements from the two 
GNSS receivers do not form a straight line, while the fusion 
of these measurements approximately makes up a straight line. 
This proves that the tested algorithm incorporating the signals 
from GNSS, and GNSS, receivers, performs correctly. 
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Fig. 4. Ship s positions measured by GNSS, and GNSS, receivers and the measurement fusion — situation 1 


Fig. 6, in turn, illustrates a situation when measurement 
data from the GNSS, receiver (more accurate) form a straight 
line, but there occurs a period when the signal fades away. The 
fusion signal, however, is maintained thanks to signals from the 
other, GNSS, receiver. This is important from the viewpoint 
of the navigational decision support system, as fading may 
destabilize the system. 

Fig. 7 presents a situation when the ship was not moving, 
lying at anchor. Also in this case for the navigational decision 
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Fig. 5. Ship s positions measured by GNSS, and GNSS, receivers and the measurement fusion — situation 2 
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Fig. 6. Ship s positions measured by GNSS, and GNSS, receivers and the measurement fusion — situation 3 
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Fig. 7. Ship's positions measured by GNSS, and GNSS, receivers and the measurement fusion — situation 4 


support system the measurement data fusion is the most proper, 
i.e. most stabilized. 

The presented situations are typical for the overall test 
results. No case was observed when the fusion of measurement 
data from a rectilinear track section would be ‘worse’ that 
autonomous signals received from GNSS receivers. 
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SUMMARY 


In this article the problem of navigational data fusion was 
analyzed. The discussed computing algorithm is based ona multi- 
sensor Kalman filter. The algorithm combines navigational data 
measurements received from single autonomous receivers. 


Implemented in the navigational decision support system, it was 
verified on the m/s Nawigator XXI. The field tests have indicated 
that the used algorithm improves the accuracy and reliability 
of measurements of navigational parameters. This confirms the 
effectiveness of algorithm operation in real conditions. 

It should be emphasized that the presented method can 
be also used in other systems, e.g. automatic control systems. 
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The Ship Handling Research and Training Centre at Ilawa is owned by the Foundation for Safety of Navigation 
and Environment Protection, which is a joint venture between the Gdynia Maritime University, the Gdansk University 
of Technology and the City of Ilawa. 


Two main fields of activity of the Foundation are: 


Ə Training on ship handling. Since 1980 more than 2500 ship masters and pilots from 35 countries were trained at 
Itawa Centre. The Foundation for Safety of Navigation and Environment Protection, being non-profit organisation 
is reinvesting all spare funds in new facilities and each year to the existing facilities new models and new training 
areas were added. Existing training models each year are also modernised, that's why at present the Centre represents 
a modern facility perfectly capable to perform training on ship handling of shipmasters, pilots and tug masters. 


Research on ship's manoeuvrability. Many experimental and theoretical research programmes covering different 
problems of manoeuvrability (including human effect, harbour and waterway design) are successfully realised at 
the Centre. 

The Foundation possesses ISO 9001 quality certificate. 


Why training on ship handling? 


The safe handling of ships depends on many factors - on ship's manoeuvring characteristics, human factor (operator 
experience and skill, his behaviour in stressed situation, etc.), actual environmental conditions, and degree of water 
area restriction. 


Results of analysis of CRG (collisions, rammings and groundings) casualties show that in one third of all the 
human error is involved, and the same amount of CRG casualties is attributed to the poor controllability of ships. 
Training on ship handling is largely recommended by IMO as one of the most effective method for improving the 
safety at sea. The goal of the above training is to gain theoretical and practical knowledge on ship handling in a wide 
number of different situations met in practice at sea. 


For further information please contact: 
The Foundation for Safety of Navigation and Environment Protection 


Head office: Ship Handling Centre: 
36, Chrzanowskiego Street 14-200 ILAWA-KAMIONKA, POLAND 
80-278 GDANSK, POLAND tel./fax: +48 (0) 89 648 74 90 
tel./fax: +48 (0) 58 341 59 19 e-mail: office@ilawashiphandling.com.pl 

e-mail: office@portilawa.com 
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